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1.  INTRODUCTION 


which  T  “efliods 

Which  can  be  used  to  dete^hie  die  se„s.,,vi„  o,  uuaier.ca.  solutions  of  sets  of 

coupl^  «„erential  equations  to  uncertainties  in  ntodeling  paranieters  entering 

prircT 

rejects  Agency  (ARPA)  under  ONR  contract  N00014-71-C-0347.  lb,s  report 
describes  work  performed  in  this  effort. 

ARPA's  interest  in  this  problem  has  been  stimulated  primarily  by  a 

need  on  the  part  of  the  defense  community  to  develop  computational  models  to 

.raor  die  certain  weapons  effects  in  which  non-eguaibrinm  chemistry  plays  an 

—  role.  Per  such  phenomena  an  accurate  treatment  of  the  hinetil  of 
the  chemical  system  requires  that: 

0)  AU  of  the  relevant  species  and  the  important  reactions 
which  occur  among  them  are  included. 

(2)  Wate  values  of  the  important  rate  coefficients  are 

For  many  applications  there  is  much  debate,  largely  based  on  intuition,  about 
i^^ich  species  are  relevant,  which  reactions  are  important,  and  which  rate  coef- 
.ente  must  be  icnown  accurately.  It  is  desirable,  therefore,  for  ARPA  m  have 

It  disposal  objective  tools  which  can  be  used  to  help  answer  guestions  of 
priority  for  the  allocation  of  their  rcsouroPQ  • 

rate  constants.  experimental  determinations  of 

A  second  motivation  tor  the  research  is  to  provide  tools  whereby  the 

stir"  “  •“'Fe  complex  systems  h,  under- 

more  u  ly  the  principles  by  which  their  constituent  parts  interact.  In  this 
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respect,  our  methods  are  not  limited  to  chemical  systems,  but  may  be  applied  to 
a  wide  variety  of  systems  for  which  computational  modeling  is  appropriate.  For 
definiteness,  however,  t^e  methods  will  be  described  in  the  context  of  non- 
equilibrium  chemical  systems  to  which  they  have  been  exclusively  applied  to  date. 

The  methods  \\liich  are  being  developed  make  use  of  the  capability  to 
solve  the  chemical  kinetics  problem  several  times  for  different  values  of  the 
rate  coefficients  within  their  expected  bands  of  uncertainty,  and  to  interpret  the 
results  of  these  calculations  in  such  a  way  as  to  determine  the  relative  sensitivity 
of  the  results  to  the  uncertainties  in  rate  coefficients.  It  is  useful  to  introduce  the 
concept  of  an  "output  function"  which  in  this  report  will  mean  a  number  which  is  of 
interest  calculated  from  the  model.  Examples  of  output  functions  are: 

•  One  of  the  concentrations  at  a  particular  time. 

•  A  concentration  maximum  over  all  time. 

•  The  time  at  which  a  concentration  reaches  some 
threshold. 

•  The  maximum  temperature  attained  by  the  system. 

•  The  radiation  intensity  at  a  particular  point  in  space 
at  a  specified  time. 

An  output  function  is  regarded  as  a  function  only  of  the  rate  coefficients  of  the 
model,  and  therefore  uncertainties  in  the  rate  coefficients  contribute  to  an  uncer¬ 
tainty  in  the  output  function. 

Since  the  dependence  of  an  arbitrary  output  function  will  be,  in  general, 
a  complicated,  non-linear  function  of  the  rate  coefficients,  traditional  analytical 
investigations  are  usually  impossible.  The  importance  of  the  various  rate  coeffi¬ 
cients  may  not  be  determined  by  investigating  each  one  separately  holding  the 
others  fixed  because  of  the  complex  interaction  between  reactions.  Furthermore, 
the  output  function  may  be  very  sensitive  to  the  rate  coefficients  in  a  localized 
region  within  the  uncertainty  limits  rather  than  at  the  upper  and  lower  bounds  of 
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uncertainty.  Such  regions  would  be  missed  by  a  sensitivity  analysis  which  examined 
only  maximum  and  minimum  values  of  the  rate  coefficients. 

Mathematically,  the  problem  reduces  to  that  of  investigating  a  non-linear 
function  of  many  variables  over  a  domain  of  its  arguments.  The  procedure  is 
therefore  characterized  by  two  parts: 

(1)  Selection  of  a  number  of  points  at  which  one  will  evaluate 
the  function. 

(2)  Determination  of  the  relative  sensitivity  of  the  function  to  the 
various  arguments  on  the  basis  of  the  function  evaluations. 

The  methods  which  have  been  developed  prescribe  a  sequence  of  values 
of  rate  coefficients  at  v^ich  the  output  function  is  evaluated.  All  of  the  rate  coeffi¬ 
cients  vary  simultaneously  along  the  sequence  of  points.  Two  different  approaches 
have  been  formulated  which  we  shall  call  the  Fourier  method  and  optimization 
methods.  In  the  Fourier  method,  the  sequence  of  rate  coefficient  values  for  the 
analysis  is  prescribed  before  any  calculations  of  the  output  function  are  obtained, 
whereas  in  the  optimization  methods,  later  points  in  the  sequence  are  determined 
by  the  values  of  the  output  function  calculated  at  earlier  points.  In  the  Fourier 
method  only  values  of  the  output  function  are  computed  vvliereas  some  of  the 
optimization  methods  require  the  calculation  of  various  partial  derivatives  of 
the  output  fimction  with  respect  to  the  various  rate  coefficients. 

The  methods  are  designed  to  accommodate  output  functions  of  an  arbitrarily 
large  number  of  rate  coefficients,  although  the  cost  of  performing  the  sensitivity 
analysis  increases  rapidly  with  this  number.  To  date  the  methods  have  been 
applied  only  to  the  relatively  simple  systems  described  below.  The  systems  are 
sufficiently  simple  that  some  results  of  the  methods  can  be  qualitatively  verified 
by  intuition  or  hindsight.  The  three  systems  which  have  been  considered  in  some 
detail  are: 


(1)  Various  "simulation  functions"  which  are  merely 
analytical  functions  of  several  parameters. 

(2)  A  chemical  system  describing  the  hi^  temperature 
dissociation  of  air  (five  concentrations,  ten 
reactions). 

(3)  A  chemical  system  describing  the  combustion  of 
hydrogen  (six  concentrations,  eighteen  reactions). 

While  these  examples  are  extremely  simple  compared  to  the  systems  for 
which  the  methods  are  needed  most,  they  do  demonstrate  the  feasibility  of  the 
methods,  and  strongly  suggest  their  applicability  to  more  complex  systems. 

In  order  to  apply  the  methods  to  chemical  systems,  one  must  be  able 
to  solve  a  set  of  chemical  rate  equations.  Many  techniques  have  been  developed 
for  such  problems,  and  a  code  called  KEM  developed  at  using  the  Gear  algorithm 
was  used.  Investigating  the  accuracy  or  efficiency  of  the  integration  technique  or 
comparing  various  techniques  was  outside  the  scope  of  the  research  effort.  The 
KEM  calculation  is  essentially  regarded  as  a  subroutine  which  gives  an  output 
function  for  a  specified  set  of  rate  coefficients. 

In  Sections  2.1  through  2.4  the  formulation  of  the  Fourier  method  is 
presented,  and  in  Section  2. 5  the  results  of  the  sensitivity  analysis  using  this 
method  for  the  two  chemical  systems  described  above  is  summarized.  In 
Section  2.6  a  Monte  Carlo  technique  for  obtaining  the  same  results  is  compared 
to  the  Fourier  method. 

hi  Section  3.1  the  optimization  methods  are  described  and  in  Section  3.2, 
comparisons  of  several  optimization  methods  are  described  for  a  simulation  func¬ 
tion.  In  Section  3. 3  an  optimization  method  is  applied  to  the  air  dissociation 
chemical  system.  In  Section  4  a  summary  of  the  work  and  proposed  extensions  are 
presented. 
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2.  THE  FOURIER  METHOD 


2. 1  FORMULATION  OF  THE  METHOD 

The  Fourier  method  is  an  attempt  to  investigate  certain  characteristics 
of  a  scalar  function  of  several  variables  over  a  convex  domain  by  simply  evaluating 
the  function  at  a  finite  set  of  points  within  the  domain.  More  specifically,  the 
function  is  evaluated  at  a  set  of  points  along  a  curved  line  within  the  domain. 

While  the  formulation  of  the  method  is  independent  of  the  nature  of  the  output 
function,  and  the  variables  on  which  it  depends,  the  notation  and  phraseology  of 
this  discussion  will  reflect  the  context  of  a  specific  application  in  which  the  output 
function,  f ,  is  a  concentration  of  a  chemical  species  in  an  isothermal,  constant 
density  reacting  system,  and  the  independent  variables  ,  i  =  1,  n  are  the  rate 
constants  which  appear  in  the  kinetic  equations  which  describe  the  reacting  system. 
The  region  of  interest  in  the  k-space  is  that  spanned  by  the  estimated  uncertainties 
in  the  current  determination  of  the  rate  constants. 

The  curved  line  on  which  f  is  to  be  evaluated  is  defined  by  the  single 
parameter,  s,  by  the  specification 

u. 

k.  =  k°  e  ' 

1  1 

and 

o  . 

u.  =  u^  sm  w.  s  ,  i  =  1,  n 

The  introduction  of  the  parameter  s  provides  a  convenient  prescrip¬ 
tion  for  varying  all  of  the  rate  constants  simultaneously.  By  selecting  a  set  of 
s  values,  a  set  of  values  of  k  are  obtained  distributed  within  the  n  dimensional 
uncertainty  domain. 


(2.1) 

(2.2) 


The  modus  operand!  of  the  Fourier  method  is  to  seleet  a  set  of  distinct 
"input"  frcqueneics  {oj.}  ,  one  for  each  reaction  rate  being  investigated.  A  set 
of  values  of  s  are  also  ehosen,  and  the  output  function  f  is  evaluated  for  each 
value  (for  the  set  of  {k^}  specified  for  each  value  of  s  throu^  Equations 
(2. 1)  and  (2.2)).  f  may  thus  be  viewed  as  a  function  of  s  , 

f(s)=  f(k^(s),  k2(s), ... ,  k^(s)]  ,  (2.3) 

and  this  function  is  Fourier  analyzed  to  obtain  the  sine  and  cosine  amplitudes 
corresponding  to  various  frequencies. 

2.2  QUALITATIVE  INTERPRETATION  OF  THE  FOURIER 
AMPLITUDES,  INTERFERENCES,  AND  ALIASING 

It  is  convenient  to  assume  for  the  moment  that  all  the  frequencies  are 
incommensurate  (irrational).  It  is  impossible  to  represent  such  frequencies 
in  a  digital  computer,  but  the  exercise  sheds  light  on  the  qualitative  interpretation 
of  the  amplitudes.  Furthermore,  it  will  be  assumed  that  f  is  an  analytic  function 
of  all  the  rate  coefficients  in  the  domain  defined  by  the  uncertainties.  If  the  ampli¬ 
tudes  corresponding  to  input  frequency  w.  and  its  harmonies  are  very  small,  the 
uncertainty  in  the  i^  reaction  rate  is  unimportant  in  determining  the  uncertainty  of  f. 
If,  on  the  other  hand,  the  amplitude  of  w.  is  large,  the  contribution  from  the  uncer¬ 
tainty  of  the  itt  rate  coefficient  to  the  output  function  is  large.  In  Section  2.4 
the  magnitude  of  the  Fourier  coefficient  will  be  related  to  an  average  of  9f/9u 
over  the  u-space  characterizing  the  uncertainties. 

In  practice,  the  unique  identification  of  a  frequency  with  the  rate  coeffi¬ 
cient  it  represents  may  be  obscured  by  two  effects:  interferences  between 
frequencies,  and  aliasing.  Interferences  arise  from  the  fact  that  incommensurate 
frequencies  may  not  be  chosen  since  they  cannot  be  represented  on  a  digital 
computer.  Commensurate  frequencies  have  the  property  that  they  cannot  be 
linearly  independent  with  respeet  to  integer  coefficients.  If,  for  example. 
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0.1+20,2=  cjg 


(li.4) 

the  amplitudes  eorresponding  to  co_  will  reflect  not  only  the  effeets  of  the  uneer- 

«5 

tainty  of  k  ,  but  also  k  and  k  through  the  interferenee. 

Aliasing  arises  whenever  a  finite  number  of  points  are  ehosen  on  the 
interval  to  evaluate  the  Fourier  amplitude,  and  the  effeets  are  more  serious  when 
evenly  spaeed  points  are  ehosen.  For  example,  if  N  equally  spaeed  points  are 
ehosen, the  amplitudes  for  wall  unavoidably  inelude  the  amplitudes  of  a  eomponent 
writh  frequeney  present  in  f(s)  whieh  satisfies  the  relation 

CO,  =  mN  -  CO,,  (2. 5) 

L  M 

where  m  is  an  arbitrary  integer. 

The  eurrent  method  takes  the  frequeneies  co.  to  be  integers  and 
evaluates  the  output  funetion  at  evenly  spaeed  diserete  values  of  s  on  the  inter¬ 
val  0  <  s  <  2t  .  This  ehoice  was  motivated  by  the  goal  of  redueing  as  mueh  as 
possible  the  number  of  points  to  be  evaluated  to  obtain  the  Fourier  amplitudes  as 
shown  in  the  next  seetion. 

2.3  CHOICE  OF  FREQUENCIES  AND  THE  METHOD  OF  FOURIER 

ANALYSIS 

The  Fourier  sine  amplitude  eorresponding  to  an  arbitrary  frequeney 

CO  assoeiated  wdth  the  Lth  rate  eoeiXieient  is  defined  as 
L  ”  T 

A  =lim-;X  |  sin  co  s  f(k  (s), . . .  ,k  (s)]ds  ,  (2.6) 

L  «iwoo  T  J  LI  n 

0 

and  the  eosine  amplitude  is  defined  as  the  same  integral  with  the  sine  function 
replaced  by  a  eosine  funetion.  Since  this  integral  is  to  be  performed  numerically 
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snd  it  is  important  to  minimize  the  number  of  s  points  used  in  the  integral,  it 
sesmt,  advantageous  to  roplaoe  this  integral  over  the  infinite  interval  by  one 
over  a  finite  interval.  Tnis  may  be  done  provided  only  rational  frequenoies  are 
present  in  the  funotior.  f .  For  output  funetions  whioh  are  suffioiently  smooth, 
one  oan  insure  that  only  rational  frequenoies  appear  by  speoifying  that  the  input 
frequenoies  themselves  are  rational.  It  is  eonvenient  to  further  simplify  the 
integral  by  tlic  transformations  in  Appendix  A  whieh  show  that  Equation  (2.6) 
with  rational  frequenoies  reduoos  to  Equation  (2.7)  with  integer  frequenoies  col 
defining  the  k's  through  Equations  (2.1)  and  (2.2), 

2n 

j*  ds  sin  u)^s  f(k^(s) . k^(s)J  .  (2.7) 

0 

integer  input  frequenoies  will  be  assumed. 

2. 3. 1  Choioe  of  Noninterfering  Frequencies 

Sinoe  integer  frequenoies  are  to  bo  input,  various  sums  and  differenoes 
of  those  frequenoies  will  ooinoide.  The  frequenoies  must  be  ohosen  so  as  to 
minimize  these  interferenoes  in  some  sense. 

We  define  an  intorferenoe  of  nth  order  as  a  linear  oombination  of  n 
frequenoies  with  ooeffioients  ±1.  Thus,  co  -  2  co„  is  an  interferenoe  of  third 

X  A 

order,  wiiile  ”  Wg  -  -  co^^  is  also  an  interferenoe  of  third  order,  that  is. 

In  order  to  understand  the  meaning  of  the  interferenoe  amplitudes  in 
relation  to  the  variation  of  the  output  funotion  with  the  rate  oonstants,  it  is  oon- 
venient  to  oonsider  f  as  a  funotion  of  the  various  u.  instead  of  as  a  funotion  of 
the  k.  through  the  transformation  in  Equation  (2.1).  Assuming  that  f  is  a 
bounded,  oontinuous  funotion  of  Uj  with  partial  derivatives  of  all  orders  in  the 


In  the  remainder  of  this  seotion. 
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region  -u.  <  u.  <  u° ,  it  may  be  expanded  in  a  Taylor  series  about  the  point 
Uj  =  0.  i  =  1,  n. 


f(u 


1’ 


u  )  = 
n' 


f  + 


i=l 


E(f..U.U  +f.  U  ) 

\  ‘J  1  j  ijj  1  j  •  7 


i  (vvj 


“k  *  'iiik  “A 


+  . . . 

where,  for  example, 

f  =  i_ 

ij  9u.9u.  21 
1  J 

The  function  f  can  also  be  expanded  in  a  Fourier  series 

CO 

f(s)  =  ^  /ajj  sin  (s  +  Bjj  cos  Is)  . 

11=1  ' 


(2.8) 


(2.9) 


(2. 10) 


9 


Substituting  u.  =  u°  sin  u.s  ii>to  Equation  (2.8)  a  series  is  obtained,  the  terms 
of  which  contain  products  of  various  powers  of  trigonometric  functions,  i.e. , 


.  p  .  q 

sin  u)jS  sin 


U).s  . 
J 


(2.11) 


By  using  the  identity 


lU).  sp 
1 
c 


iu.sq 
e  J 


is(po.  +  qu  ) 

i  J 


(2. 12) 


and  realizing  that  pw. +qu)^  is  integer,  one  observes  that  the  two  series  contain 
exactly  the  same  terms  making  the  identification  of  the  Aj^,  Bjj  with  various  sums 
of  partial  derivatives  evaluated  at  the  point  u.  =  0,  i  =  1 ,  n.  Furthermore,  a  p  th 
order  interference  in  the  Fourier  series  arises  from  terms  in  the  Taylor  expansion 
of  order  p  or  greater.  In  fact,  if  f  were  linear  in  all  u's  ,  the  only  frequencies 
appearing  the  Fourier  decomposition  would  be  the  input  frequencies. 


The  assumption  is  tacitly  made  that,  in  general,  high  order  interferences 
(and,  therefore,  hi^  order  terms  in  u)  are  smaF  compared  with  lower  order  inter¬ 
ferences  and,  therefore,  it  will  be  adequate  to  select  the  input  frequencies  in  such  a 
way  as  to  eliminate  effects  of  interferences  below,  say,  order.  Interferences 
above  order  will  give  rise  to  inaccuracies  in  the  amplitudes.  From  the  stand¬ 
point  of  interferences,  therefore,  one  prefers  to  choose  frequencies  consistent  with 
large  M ,  which  turn  out  to  be  relatively  large  int^ers.  On  the  other  hand,  the 
number  of  output  fimetion  evaluations  necessary  to  evaluate  the  Fourier  amplitudes 
of  large  frequencies  is  greater  than  that  necessary  for  small  frequencies.  The 
method  which  has  been  investigated  to  date  makes  use  of  frequencies  for  which 
no  interferences  of  the  order  of  less  than  five  coincide  with  the  input  frequencies. 
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2. 3.1.1  Numerical  Piocedure  for  Determining  the  Frequencies 

The  method  currently  used  to  determine  the  integer  input  frequencies, 
u. ,  is  essentially  a  trial  and  error  procedure.  The  object  is  to  find  a  set  of 
;i  integers  co.  such  that  no  ii  terferences  of  the  order  of  four  or  less  will  coincide 
with  an  iiput  frequency.  Explicitly, 


pWj  +  qwj^  +  rwjj  +stj^?^u). , 

where  . cOj^,  co^,  and  w.  are  input  frequencies,  and  thr.  coefficients 
obey 

-4  <  p ,  q,  r,  s  <  4  , 
and 

|p|+|q|  +  lr|+|s|s4  .  ,2.13, 

An  obvious  simplification  is  made  by  chosing  all  odd  integers  for  the 
input  frequencies  so  that  there  are  no  coincidences  with  second  order  or  fourth 
order  interferences.  The  only  conditions  that  remain  to  be  satisfied  are  that 


and 


U>.  +U)  ^  u  -  (J. 


(2. 14a) 


(2. 14b) 


The  first  step  is  to  find  an  ordered  set  of  differences  such  that  Equa¬ 
tion  (2. 14a)  is  satisfied,  that  is,  a  set  of  d.  such  that 


i=l 


(2. 15) 


1  1 


are  distinct  for  1  <  m  <  n  - 1 .  This  is  done  recursively,  and  while  there  exist 
many  solutions  for  a  particular  n .  one  wishes  to  find  the  set  with  the  minimum 

largest  difference  d^^^ 


n-1 

d  =  d. 

max  tm-i  I 

i=l 


(2. 16) 


which  determines,  in  part,  the  largest  frequency. 

The  process  proceeds  as  follows.  A  table  is  begun  for  n  -  3 ,  with  the 
ordered  set  of  frequencies  <  u)2<  Wg .  and  the  ordered  set  of  differences, 
di  =  Ug  -  .  and  d^  =  iJg  “  Wg-  d.  are  determined  to  be  the  smallest 

possible  even  numbers  (excepting  2),  i.e. , 


The  difference  d*  =  d^+ d^  =  10  also  occurs  and  is  distinct  from  d^  and  d^  . 

It  is  convenient  to  express  the  results  in  a  triangular  table 

4  6 

10 

To  continue  the  process  for  n  =  4 .  one  chooses  dg  =  8  (must  be  distinct)  and  the 
d*  are  determined  by  forming  the  diagonal  as  shown  below 

i  =  1  2  3 

4  6  8 

10  14 
18 

from  upper  right  to  lower  left  (i  =  n-l.  i  =  l).  the  element  in  the  ith  column  being 


the  sum  of  the  diagonal  in  the  i  +  lst  column  and  the  difference  appearing  at  the 
top  of  the  column.  The  successive  d.  are  chosen  so  that  no  number  in  the 
ei.Mre  tabic  is  repeated. 

One  is  not  restricted  to  choosing  the  next  d.  to  be  the  lowest  unused 
even  integer,  but  experience  shows  that  as  a  general  rule  this  choice  will  result 
in  the  lower  frequencies.  Two  sets  of  differences  have  been  calculated  and  are 
presented  in  Tables  2. 1  and  2. 2.  The  set  of  differences  in  Table  2. 1  use  the 
lowest  possible  unused  difference  to  increment  n  ,  while  in  Table  2.2  d^  was 
taken  not  to  be  the  lowest  unused  even  number,  6,  but  8.  It  appears  that  the 
frequencies  from  either  method  are  equally  good,  and  that  the  results  are  not 
strongly  dependent  on  the  ordering  of  small  differences.  Furthermore,  since 
the  even  numbers  between  dj^  and  d^  which  are  omitted  is  relatively  small, 
it  appears  that  little  improvement  could  be  madf  by  finding  the  optional  ordering 

since  d  would  not  change  significantly. 

max 

The  set  of  differences  for  n  ordered  frequencies  is  the  set  of  the  first 
n  -  1  differences  of  th*^  top  row  in  either  table.  For  example,  for  n  =  6,  the 
differences  in  Table  2. 2  are:  4,  8,  6,  10,20. 

After  a  set  of  differences  are  found  for  a  particular  n ,  the  whole  set 
of  frequencies  is  detex’mined  by  specifying  the  lowest  frequency.  In  the  example 
above,  with  n  =  6  ,  one  assigns  the  differences  in  reverse  order  and  assigns 
=  1  to  obtain  the  frequenej'  set:  1,  21,  31,  37,  45,  49.  The  second  step  is  to 
add  to  all  frequencies  the  same  integer  until  Equation  (2. 14b)  is  satisfied, 

U).  +  U.  /  (0,-0)/)  .  (2* 

1  ]  k  X 

Frequencies  for  which  no  interferences  of  order  less  than  5  coincide  with 
an  input  frequency  have  been  calculated  for  n<  19  and  the  frequency  sets  are 
shown  in  Table  2.3. 
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Table  2.3 


FREQUENCY  SETS  AND  CORRESPONDING  NUMBER  OF 
POINTS  AVOIDING  INTERFERENCES  AND 
ALIASES  OF  INTERFERENCES  THROUGH  4^  ORDER 


Dimensionality  Frequency  Set  Minimum  Number  of  Points 


5 

11,  21,  27,  35,  39 

142 

5* 

2,  42,  62,  74,  90 

191 

6 

1,  21,  31,  37,  45,  49 

182 

6* 

2, -42,  62,  74,  90,  98 

231 

7 

17,  39,  59,  69,  75,  83,  87 

334 

8 

23,  55,  77,  97,  107,  113, 

121,  125 

486 

9 

19,  59,  91,  113,  133,  143 

149,  157,  161 

630 

10 

25,  63,  103,  135,  157,  177, 
187,  193,  201,  205 

806 

11 

41,  67,  105,  145,  177,  199, 
219,  229,  235,  243,  247 

974 

12 

31,  87,  113,  151,  191,  223, 
245,  265,  275,  281,  289,  293 

1158 

13 

23,  85,  141,  167,  205,  245, 
277,  299,  319,  329,  335,  343, 
347 

1374 

14 

87,  133,  195,  251,  277,  315, 
355,  387,  409,  429,  439,  445, 
453,  457. 

1814 

15 

67,  143,  189,  251,  307,  333, 
371,  411,  443,  465,  485,  495, 
501,  509,  513 

2038 

16 

73,  169,  245,  291,  353,  409, 
435,  473,  513,  545,  567,  587, 
597,  603,  611,  615 

2446 

17 

85,  145,  241,  317,  363,  425, 
481,  507,  545,  585,  617,  639, 
659,  669,  675,  683  ,  687 

2734 

18 

143,  229,  289,  385,  461,  507, 
569,  625,  651,  689,  729,  761, 
783,  803,  813,  819,  827,  831 

3310 

19 

149,  275,  361,  421,  517,  593, 
639,  701,  757,  783,  821,  861, 
893,  915,  935,  945,  951,  959, 

3848 

963 

^Tliese  sets  of  S  and  6  frequencies  are  not  minimal,  but  have  been  used  and  are 
only  included  for  completeness.  They  are  also  free  of  Inteferences  and  aliases 
of  order  less  than  S. 


2.3.2  Fourier  Analysis  and  Avoidance  of  Aliases 


The  Fourier  analysis  could  be  performed  either  by  evaluating  f  at  N 
evenly  spaced  or  unevenly  spaced  points  on  the  interval  (0,  2ir).  Evenly  spaced  points 
are  chosen  although  it  is  not  certain  that  it  is  an  optional  choice.  Given  evenly 
spaced  points,  howevc- ,  it  is  convenient  to  replace  the  integral  in  Equation  (2.7) 
by  a  sum.  Consider  the  sum 

E  '(^)  •  <2-“) 

i=i 

It  is  shovwi  in  Appendix  B  that 


A*  =  A  +  error  , 
Jj  Li 


(2. 19) 
/ 


where  the  error  terms  called  "alias  amplitudes"  are  the  amplitudes  of  various 
interference  frequencies: 


00 

error  = 

m=0  j 


(mN-  a\) 


where 


u. 

1 


is  any  frequency  appearing  in  f  satisfying  the  relation 


(2. 20) 


mN-  w.  =  CO  I  (2.21) 

1  ^ 

where  m  is  an  arbitrary  integer.  The  number  N  is  chosen  so  that  Equation  (2.21) 

is  not  satisfied  for  u  an  interference  of  order  less  than  5  and  w  is  any  input 

J  3j 

frequency.  This  number  is  determined  by  calculating  for  a  trial  N" 

mN  -  co_ 
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where  is  an  input  frequency,  and  comparing  this  number  with  each  inter¬ 
ference  frequency  through  fourth  order.  N  is  the  minimum  value  of  N  for 
which  there  are  no  coincidences.  Tlie  number  N  for  the  frequency  sets  in 
Table  2. 3  is  also  shown. 


2.4 


QUANTITATIVE  INTERPRETATION  OF  THE  FOURIER 
ALTERNATE  WEIGHT  FUNCTIONS  IN  u  SPACE 


AMPLITUDES; 


One  can  view  the  Fourier  method  as  being  characterized  by  a  sampling 
prescription  in  n-dimensional  u  space  in  addition  to  an  anali^sis  of  the  output 
function  evaluated  at  the  selected  points.  The  result  of  the  method  depends 
of  all  the  points  sampled,  and  we  show  below  that  certain  Fourier  amplitudes 
are  related  to  averages  of  the  output  function  with  various  weight  functions  over 
the  n-dimensional  u  space. 


It  is  shown  in  Appendix  C  that  if  f  is  a  polynomial  funption  of  the  u’s 
of  order  p-1 ,  and  the  input  integer  frequency  set  is  chosen  so  that  no  interferences 
of  order  less  than  p  coincide  with  input  frequencies,  then  the  one-dimensional 

integral  corresponding  to  the  Sine  amplitude  of  input  frequency  is  equal  to  an 

n-diniensional  integral  over  u  space: 


(This  equality  does  not  hold  for  functions  f  of  order  higher  than  p-1  and  an 
error  term  appears  proportional  to  (f^V2P'^  as  shown  in  Appendix  C. )  This 
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Fourier  amplitude  is  therefore  equal  to  the  L ^  eomponent  of  the  eentroid, 
weighted  by  the  funetion 


An  altei'native  form  of  the  multiple  integral  may  be  derived  by  integrating  by 

parts  with  respeet  to  u  ,  i.e. , 

L 


The  weight  fimetions  in  these  multiple  integrals  emphasize  various  regions 
of  u-space.  In  the  form  in  Equation  (2. 22),  all  edges  of  the  hypereube  are 
emphasized  equally,  whereas  in  the  integral  in  Equation  (2.24),  the  distribution 
is  uniform  in  u^^  but  emphasizes  the  edges  in  the  other  dimensions.  Efforts 
to  modify  this  weighting  will  be  deseribed  below. 

The  implieation  of  Equation  (2. 24)  is  that  the  Fourier  eoeffieient  may  be 
interpreted  as  an  average  over  the  Anijole  uneertainty  region  of  9f/9  u^  ,  and  it 
therefore  represents  a  measure  of  sensitivity.  Alternatively,  by  analyzing  inf 
instead  of  f  one  obtains  the  measure  of  sensitivity  of  the  relative  uneertainty  in 

f  to  the  uneertainty  of  the  various  rate  eoeffieients,  i.  e. ,  an  average  of  -  — 

f  9u  * 
L 

The  estimated  experimental  uneertainty,  u°,  represents  the  half  width 
of  some  distribution  of  probable  values  for  the  reaetion  rate.  Therefore,  one 
ean  argue  that  the  sampling  procedure  should  emphasize  the  region  of  u-space 
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o 


o 


u 


) 


near  the  origin.  The  appearanee  of  the  wei^t  functions  which  emphasize  the 
boundaries  of  the  region  may  therefore  be  undesirable.  One  can  alter  the 
weight  function  by  choosing  an  alternative  definition  of  u.  (sec  Equation  (2.2): 


u.  =  u?  V.  (sin  u.  s  J  ,  (r. 

1  1  1  i  i' 

and  choosing  v.  to  be  a  monotonie  function  of  sin  o)  s  . 

‘  i 

We  will  demonstrate  this  by  finding  Vj  such  that  the  weighting  function 
in  Equation  (2.22)  is  given  by 


n 


W(u)  = 


n 

i=l 


2  -  o2 


(2.26) 


Assuming  that  f  is  a  polynomial  in  the  u's  of  order  p-1 ,  and  using  the 
same  argument  as  those  in  Appendix  C,  one  ean  show  that 


1 

(27r) 


.  (2.27) 


Setting 


du. 

1 


_1_ 

o 

u. 

1 


ds, 


dv. 


GC 


2/  o" 

1 

“T  e 


(2. 28) 


) 
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one  may  solve  for  v.(s.)  .  A  solution  for  this  equation  for  -1  <  v.  <  1  is 
v.(b.)  =  erf’^  sin'^  b  J  • 


(2.29) 


where 

b.  =sina;.s.  (2.30) 

and 

-tt/2  <  sin  ^  bj  <  V2  . 


This  gives  for  the  right  side  of  Equation  (2.27) 

^ r\  r\  %t  ^ 


(2erf(l)) 


o  '<“r“2 . V- 


(2.31) 


-u. 


-u 


L  i=l 


One  ean  integrate  this  by  parts  over  u  to  obtain  an  average  of  9f/9u_  with  a 

Ju  L 

Gaussian  weight  funetion  remaining  over  all  u’s  exeept  u  ,  and  a  weight 

L 

funetion  in  the  u  dimension  proportional  to 
Ju 

H(Ul)=  Jdu^^u^e 
-00 

whieh  is  also  peaked  near  u  =  0 . 

Li 

Some  ealeulations  reported  in  the  next  seetion  used  this  formulation  of  the 
Fourier  method,  but  the  results  do  not  differ  significantly  from  those  obtained 
with  the  original  formulation. 


u. 


/ 


du’  u’  e 

Ju  Li 


(2.  32) 


21 


o 


o 


2. 5  APPLICATIONS  OF  THE  FOURIER  METHOD 


O 


o 


o 


o 


o 


The  Fourier  method  has  been  applied  to  various  simulation  functions 
(analytical  functions)  and  to  results  calculated  from  models  of  two  chemical 
systems.  The  results  for  the  chemical  systems  will  be  presented  in  this 
section. 

2.5.1  The  Nitrogen-Oxygen  (N-0)  System 

The  N-0  system  is  modeled  by  the  following  set  of  chemical  reactions, 
rate  constants  and  equilibrium  constants: 


Og  +  M  ^20  +  M 

N2  +  M522N  +  M 

NO  +  M  ^N  +  0  +  M 

Ng  +  0  JSNO  +  N 

^4. 

0  +  N  5=  NO  +  0 

and  is  characterized  by  the  following  rate  equations, 


0 


dCf 

^  =  e^M)  .k^(c^c^.^  c^) 


D 


D 


'5C2) 

(cjM-^c^c^m)  -k^(c302 


■k,  "^5 '4)- 
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(It  "  "*^2  (°3  ^  ■  K  °4  (°3  °2  "  K  °5 


=  21^2  (cgM-^  c^m)  ^  c^c^m) 


■^*^4  Iv2-k:°4 

'  4 


) ->'E  (v4-ir  V2)  • 


=  -  '‘3  (“5“  -  #3  V2“)  '‘4  (‘=3“2  ■  '5°4) 

*“5  (V4  -ib  V2) 


(2.33) 


where 


“l  =  I°2l 

C3  =  (0) 

<=3  =  '“2> 
»4  =  (N) 

O3  =  (NO) 


=  E 


(2.34) 


The  assumption  of  equal  reaction  rates  with  each  species  acting  as  a  third  body, 
M,  for  reactions  1,  2,  and  3  simplifies  the  system,  k.  are  the  rate  constants 
for  the  reaction  in  the  forward  direction  as  written,  and  K.  are  the  equilibrium 
constants.  The  rate  constants  for  the  reverse  reactions  are  related  to  k.  by  the 
relation 
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i  reverse 


(2.35) 


''i  =  k 


The  system  was  assumed  to  react  at  constant  temperature,  6000°K,  and 

constant  density  beginning  at  t  =  0  with  only  N  and  O  present.  The  reaction 

2 

rates  and  equilibrium  constants  used  and  the  initial  conditions  assumed  are 
shown  in  Table  2. 4.  A  plot  oi  the  time  evolution  of  this  system  is  shown  in 
Figure  2.1. 

This  chemical  system  is  interesting  because  it  has  in  addition  to  the 

monotonic  concentrations,  the  NO  concentration  which  passes  through  a  maxi- 
-4 

mum  near  t  =  10  sec. 


The  rate  constants  in  Table  2.4  represent  experimentally  obtained  rate 
coefficients,  and  while  the  uncertainties  of  these  particular  rate  constants  are 
not  large,  it  was  arbitrarily  assumed  that  the  experimental  uncertainties  for 
all  the  coefficients  was  ±2  orders  of  magnitude.  In  the  sensitivity  analysis, 
values  of  the  rate  coefficients  were  investigated  lying  in  the  range 


k 


i 


±4.606 

e 


i  =  1,  5 


(2. 36) 


Throughout  the  course  of  the  study,  various  sets  of  frequencies  were 
used,  some  were  not  the  minimal  set  described  earlier,  but  all  sets  had  the 
property  that  no  interferences  of  the  order  less  than  5  coincide  with  a  fundamental. 
The  most  common  set  of  frequencies  used  was 

(2,  42,  62,  74,  90) 

for  vhich  the  lowest  number  of  evenly  spaced  points  viiich  avoids  aliasing  of 
fourth  order  interferences  is  191. 


Table  2.4 


PARAMETERS  AND  INITIAL  CONDITIONS 
FOR  CALCULATIONS  OF  THE  N-0  SYSTEM 

Initial  Conditions 
[N  1  =  8  X  10  moles/cc 

A 

[Oj  =  2  X  10"®  " 

Nominal  Values  of  the  Rate  Coefficients 
k°  =  8. 5  X  10^®  (moles/cc)"^  sec"^ 

k°  =  3.0xlo'^  " 

kg  =  8.0x10®  " 

k°  =  9.0xl0^®  " 

4 

k°  =  8.0xl0^^  " 

5 

Equilibrium  Constants 
=  7.8  X  lO"^  (moles/cc) 

Kg  =1.0x10"^ 

Kg  =  1.3x10"®  " 

K  =9.0x10"® 

4 

K-=  5.9x10^  -  '  ^ 

_ _  L 
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Parameterizing  the  rate  constant  parameters,  u  by  the  relation 


u.  =  4.606  sin  w.  s 
1  1 


(2. 37) 


and  varying  s  in  the  range  0  <  s  <  2Tr,  the  concentrations  of  the  various 
species  were  calculated  by  integrating  the  rate  Equations  (2.2)  to  a  particular 
time  for  each  of  the  191  s  values.  The  191  vslues  of  a  concentration  at  a  partic¬ 
ular  time  defines  an  output  function  which  may  be  Fourier  analyzed  by  computing 
the  Fourier  coefficients  corresponding  to  the  input  frequencies  according  to  formula 


N 


A  -1  •  277  jgL  ,/27ri\  M-101 

N  2-i  N  N  /  ’  ^ 


(2. 38) 


f  is  the  output  function  evaluated  at  s  =  .  ITie  relative 

magnitudes  of  these  amplitudes  is  taken  to  be  a  measure  of  the  rielative 
sensitivity  of  the  uncertainty  in  the  output  ftmction  to  uncertainty  in  the  cor¬ 
responding  rate  coefficients. 


As  an  example  of  f(s),  Figure  2.2  shows  a  plot  of  the  NO  concentration 

-4 

as  a  function  of  s  at  t  =  10  sec.  The  curve  consists  of  points  corresponding 
to  the  NO  concentration  for  each  discrete  value  of  s ,  v\diich  have  been  connected 
by  straight  lines.  The  Fourier  coefficients  of  this  function  are  presented  in 
Section  2.5. 1.2. 


2. 5.1.1  Interpretation  of  Results  for  the  N-0  System 
at  10“6  sec 

mmm*  A 

The  N-0  system  was  run  to  a  time  of  lO”  sec  and  the  Fourier  analysis 
done  on  the  concentrations  at  that  time.  The  results  for  this  calculation  are 
presented  in  Table  2.  5. 
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Table  2. 5 

SENSITIVITY  OF  THE  CONCENTRATIONS  AT  lO"®  SEC 
TO  UNCERTAINTIES  IN  THE  RATE  COEFFICIENTS  FOR 
THE  N-O  SYSTEM 


Input  Frequency 

Corresponding  Rale  Coefficient 

Sensitivity 

(Ojjl  X  10® 

2 

’'l 

-4.47 

90 

‘e 

-0.0244 

V4 

-0.0227 

42 

-0. 00047 

62 

-0. 00039 

-  (oj  X  10®  -  -  -  . 

2 

•'l 

8.53 

74 

‘4 

-0. 36 

90 

0. 018 

62 

-0.  029 

42 

0.0027 

-  (Njj)  X 10®  - 

2 

‘'l 

-0. 389 

74 

-4 

-0. 386 

42 

‘a 

-0.  0087 

62 

-0.  0039 

90 

-0. 00103 

- •  (N)  X 10®  - - 

2 

‘'l 

0.367 

74 

-4 

0.362 

90 

-0.  028 

62 

-0. 013 

42 

‘a 

-0. 011 

-  (NOj  xlO^®  .  - 

2 

>^1 

0.412 

74 

*'4 

0.410 

90 

*'6 

0.031 

42 

-0. 0066 

62 

•'3 

-0.0048 

29 


In  (he  column  labeled  "SenelOvlty"  are  listed  the  Fourier  amplitudes 
corresponding  to  (he  rate  coefficient  shown.  At  this  early  time  (see  Figure  2. 1) 
onfy  reaction  1  is  important  in  the  production  of  O  from  and,  (hus,  only 
the  uncertainty  m  contributes  to  these  concentrations.  The  sign  of  the 
amplitude  may  be  understood  by  Equation  (2.24)  showing  the  amplitude  proportional 
to  an  average  of  8f/au.  Thus,  the  negative  sign  in  the  sensitivity  of  O  to  k 
means  that  8f/8kj  is  negative,  as  one  would  expect.  Since  k  is  small  the  ' 
primary  source  of  N  is  through  reaction  4.  One  can  show  that  at  these  early 
times, when  the  and  are  not  substantially  depleted,  the  concentrations 
of  NO  and  N  wUl  be  proportional  to  the  product  k^k^ .  The  Fourier  coefti- 

eients  of  these  species  are  therefore  (approximately)  symmetric  in  k  .k  . 

1  4 

Ii  is  desirable  that  the  results  are  jidependent  of  the  input  frequencies. 
Small  amplitudes  can  depend  strongly  on  the  input  frequency  set  through  inter¬ 
ferences  obscuring  the  fundamentals.  This  provides  a  means  for  estimating  the 
effects  of  interferences  and  aliases,  for  example,  one  can  reassign  the  frequencies 
to  the  various  rate  constants  and  repeat  the  analysis.  If  the  relative  sensitivity 
is  invariant  under  permutations,  the  result  is  reliable.  If  it  varies  with  permuta¬ 
tion,  the  result  is  questionable.  For  the  amplitudes  in  Table  2. 5  only  the  smallest 
amplitudes  will  be  unreliable.  A  comparison  of  results  under  permutation  is  pre- 
sented  below  for  another  set  of  output  functions. 

2.5.1.2  ^“‘.YPretntion  of  Besults  for  the  N-0  System  at 

The  N-0  system  with  the  same  initial  conditions  and  rate  constant  speci¬ 
fication  was  run  to  t  =  10-  sec,  and  the  various  concentrations  were  analyzed 
With  the  results  shown  in  Table  2. 6. 

Since  the  time  is  later,  more  reactions  (and  more  uncertainties)  have 
come  into  play.  If  we  arbitrarUy  choose  a  factor  of  10  to  separate  the  "important" 
from  (he  "unimportant"  rate  uncertainties.  [Oj]  and  [o]  are  sensitive  only  to 
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Table  2.6 


SENSITIVITY  OF  THE  CONCENTRATION  AT  lO"^  SEC 
TO  UNCERTAINTIES  IN  THE  RATE  COEFFICIENTS 
FOR  THE  N-0  SYSTEM 


Input  Frequency 


Corresponding  Rate  Coefficient 


Sensitivity 


tOjj]  X  10® 


t0)x  10 


-0.00146 


iNg)  X  109 


[N]  X  10* 


[NO)x  10* 
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the  uncertainty  in  .  On  the  other  hand,  the  other  concentrations  depend  on 
several  uncertainties.  It  is  interesting  that  the  most  direct  route  to  produce 
N ,  reaction  2,  is  relatively  unimportant  in  determining  any  concentration. 


In  order  to  investigate  possible  errors  in  these  results  due  to  frequency 
interference  effects,  an  additional  analysis  was  made  changing  the  frequencies 
assigned  to  the  various  rate  coefficients.  A  cyclic  permutation  of  the  frequencies 
was  used  which  caused  different  interferences  of  order  >  5  to  coincide 
with  the  fundamentals.  The  problem  for  which  this  comparison  was  made  was  the 

same  as  that  described  in  Table  2.4  except  that  one  equilibrium  constant,  K  , 

-1  ° 
was  changed  to  5. 9  x  10  which  changed  the  reverse  reaction  rate,  k_ .  The  re- 

5 

suits  of  this  calculation  are  presented  in  column  II  of  Table  2.7.  The  problem 
with  modified  Kg  was  also  run  using  the  same  frequency  set  as  that  used  in 
Table  2.6  and  the  results  presented  in  the  column  labled  "I"  so  that  direct  compar¬ 
ison  is  possible.  Only  the  small  amplitudes  are  changed  significantly.  In  no  case 
did  the  sign  change  due  to  permutation  although  the  rank  based  on  the  magnitude  of 
the  various  small  amplitudes  changes  occasionally. 

In  Section  2.4  an  alternative  parameterization  of  k^  in  terms  of  the 
parameters  was  discussed  which  resulted  in  a  Gaussian  weight  function  for  the 
multi-dimensional  integral.  In  columns  III  and  IV  the  results  for  the  same  N-0 
system  with  modified  Kg  are  presented,  with  the  original  frequency  set  used  in 
colunm  in  and  the  permuted  set  used  in  column  IV.  Again,  significant  differences 
due  to  permutation  occur  only  in  the  relatively  small  coefficients.  It  is  interesting 
to  compare  the  results  in  columns  I  and  in  due  to  the  different  sampling  of 
points.  The  formulation  in  column  I  favors  points  near  the  extremes  of  the 
uncertainty  ranges.  The  general  tendency  for  the  magnitude  of  the  amplitudes  to 
be  slightly  smaller  in  column  III  indicates  that  af/3u  tends  to  increase  with  u 
in  general.  There  do  not  seem  to  be  significant  differences  in  the  results  between 
parameterizations  of  k^  .  Moreover,  the  more  complex  parameterization  will,  in 


Table  2.7 


EFFECTS  OF  PERMUTATION  OF  FREQUENCIES  AND  ALTERNATE 
WEIGHT  FUNCTIONS  ON  SENSITIVITY  ANALYSIS 
CONCENTRATIONS  AT  SEC  FOR  THE  N-0  SYSTEM 


Rite  Coefficient 

1 

U 

ID 

IV 

. '°2» 

X  10*  ... 

*1 

-1.14 

-1.14 

-1.05 

-1.04 

*6 

-0.040 

-0.038 

-0.0098 

-0.036 

^3 

-0.040 

-0.015 

-0.0248 

-0.047 

^2 

-0.028 

-0.038 

-0.0097 

-0.016 

*4 

-0.018 

-0.020 

-0.0030 

-6.030 

. |OI  X 

10*  .... 

V 
•  1 

2.27 

2.28 

2.08 

2.05 

0.124 

0.127 

0.074 

0.111 

'‘s 

0.105 

0.097 

0.031 

0.075 

^2 

0.045 

0.069 

0.015 

0.034 

-0.011 

-0.039 

-0,074 

-  0. 152 

. iNj)  * 

10* . 

•0.569 

-0.555 

-0.428 

-0,363 

,  ^ 

•0,465 

-0.456 

-0.292 

-0.246 

•0.290 

-0.309 

-0.159 

-0.146 

•0,069 

-0.0<!6 

-0,035 

-0.058 

^2 

•0.045 

-0,078 

-0.023 

-0.022 

1.12 

1.09 

0.836 

0.693 

0.653 

0.838 

0.504 

.  0,400 

0.597 

0.640 

0.360 

0.331 

*'3  • 

0.171 

0.169 

0.093 

0.137 

0.079 

0.149 

0.012 

0.016 

7.71 

7,76 

7.99 

9.28 

^3 

•3.31 

-3.70 

2.40 

2.37 

^5 

•1.71 

-2.17 

4.07 

3.96 

1.40 

1.58 

2.00 

3.33 

kjj 

1.06 

0.721 

0.395 

0.316 

I  Soiiltlvltjr  with  frc(|ucnc7  ict  <2.  42.  62,  74,  ?0)  rorrcspondine  to  (kj . k 

H  ^litivity  with  permuted. frequencies  (62,  74,  90,  2,  42)  correnpondine  to 

’’‘l . V’ 

m  Sensitivity  using  sllornatc  form,  Equstlons  (2, 25,  snd  2.29)  with  frequency  set 
(2.  42  ,  62.  74  .  90).  . 

^t  ltlvlty  usini;  siternste  form,  Equstlons  (2,25  snd  2.29)  with  frequency  set 
(62,  74  ,  90.  2,  42), 


IV 


s”'  ‘  wIT  “r  «  ‘“O  "Ot  to  for  a.e  N-0 

tem).  We,  therefore,  tend  to  favor  the  staple  preempt, on. 

The  reenlte  of  the  eeneltivlty  analysis  -or  this  simple  N-0  system 
3gree  with  intuition  in  the  cnsps  urViav>v>  •  *,  -a* 

_  .  .  ®  ‘htuition  gives  an  answer.  Furthermore 

ambiguity  In  the  results  due  to  freauenev  Int.rf  raermore, 

1  .  .  irequency  interferences  Is  small  for  determining 

tntlC  ■*“  “““  «nsi,ivit,es  may  he  quant,- 


Oil 


2.5.2  The  Hydrogen-Oxygen  System 


The  H-0  system  is  modeled  by  the  following  set  of  chemical  reactions: 


H  +  2  H  +  M 

2 

•"i 

0+M«20+M 

2 

^2* 

'^2 

OH  +  0  +  H  +  M 

HO  +  M3£H  +  oh  +  M 

2 

H  +  0  ^OH  +  0 

2 

^5’ 

0  +  H  OH  +  H 

2 

^6* 

H  +  OH5SH  0  +  H 

2  2 

^7 

00 

OH  +  OH»H  0  +  O 

^9’ 

The  system  was  assumed  to  react  at  constant  temperature,  2000%  and 
constant  density  beginning  with  only  H  and  O  present.  The  reaction 

A  2 

rates  and  equilibrium  constants  used  and  the  initial  conditions  assumed  are  shown 
in  Table  2. 8.  Both  forward  and  backward  reactions  were  considered  so  that 
there  are  9  independent  uncertainties. 

The  early  and  late-time  concentrations  are  shown  in  Figures  2. 3  and 
2.4.  In  Figure  2. 3  the  concentrations  of  H  and  O  \(4iich  are  essentially 

M  2 

unchaiged  from  the  initial  conditions  have  not  been  shown.  At  early  times  the 

only  reactions  v^ich  contribute  are  numbers  1,  2,  and  8  since  no  atomic  O 

or  H  or  OH  is  initially  present.  As  more  O,  H,  and  OH  become  available, 

reactions  5,  6,  and  7  become  dominant  which  results  in  the  sharp  decrease 

in  o^grgen  and  hydrogen  molecules  and  the  correspondingly  sharp  increase  of 

"3 

other  species  at  t~10  sec.  Finally,  the  slower  reactions  equilibrate 
slowly  until  a  relative  steady  state  is  reached  at  about  t  =  lo"^  sec.  This 
steady  state  is  not  true  equilibrium  in  the  sense  that  reactions  2,  5,  and  8 
are  still  quite  far  from  equilibrium,  but  further  relaxation  to  equilibrium  is 
extremely  slow. 
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Table  2. 8 


PARAMETERS  AND  INITIAL  CONDITIONS  FOR 
CALCULATIONS  OF  THE  H-0  SYSTEM 

Initial  Conditions; 

(1^2^  “  0  X  10  moles/cc  [O^]  =  4, 0  x  10  ^  moles/cc 

Nominal  Values  of  the  Rate  Coefficients; 

=  5. 78  X  10  (moles/cc)  ^  sec  ^  =  4.0  x  10^^  (moles/cc)~^  sec”^ 


k°  =  4.47x10^ 

A 

ft 

1  0  ,  ,  «13 

ky  =  1. 51  X  10 

ff 

kg  =  1. 03  X  10® 

If 

kg  =  2.20x10® 

ft 

k°  =  6.79xl0® 

4 

ft 

k°  =  1.6x10^^ 

9 

ff 

k°  =  1.23x10^^ 

5 

ft 

Equilibrium  Constants: 

=  1.59x10”^^ 

moles/cc 

Ko  = 1.4  - 

Kg  =  2. 69  X  10"^^ 

If 

K^  =  9.76  - 

Kg  =  1. 14  X  10"^^ 

If 

Kg  =  2. 02  - 

K^  =  1.63X  10“^^ 

If 

Kg  =  6. 76  “ 

K^  =  2.37  X  10“^ 

0 

- 
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The  rate  constants  shown  in  Table  2. 8  are  taken  to  be  the  'Tbest  guess” 
values  and  one  order  of  magnitude  uncertainty  in  each  rate  constant  was 
assumed  in  order  to  demonstrate  the  technique.  The  frequencies  for  the 
analysis  were  taken  to  be  the  minimal  set, 

w  =  (19,  59,  91,  113,  133,  143,  149,  157,  161) 

630  evenly  spaced  points  were  taken  on  the  s  interval 

0<s<  2jr 

to  evaluate  the  Fourier  amplitudes.  The  aualysis  was  run  on  the  concentrations 
at  one  early  time  and  one  late  time. 

2. 5. 2.1  Interpretation  of  Results  for  the  H-O  System  at  10~^  sec 

The  results  for  the  sensitivity  analysis  for  t  =  lo"®  sec  are  shown  in 
Table  2. 9.  At  this  very  early  time,  the  and  concentrations  have  not 
varied  from  their  initial  values  so  that  the  sensitivity  analysis  showed  their 
Fourier  amplitudes  to  be  very  small,  due  only  to  roundoff  in  the  analysis. 

The  relative  sensitivites  of  the  other  species,  however,  reflect  the  importance 
of  the  various  reactions  and  the  uncertainties  in  their  rate  constants.  In  the 
case  of  H,  reaction  1  is  the  prime  contributor  and,  therefore,  its  rate 
coefficient  is  crucial  in  determining  the  accuracy  of  [H]  .  k  and  k  are  also 

important  since  reactions  7  and  8  together  produce  H.  ThelensitivLs  of 

the  two  reactions  are  the  same  since,  at  early  times  these  two  reactions  yield 

[H]'=^k^kg  , 

and  since  the  uncertainties  u°  =  u°  ,  the  results  are  symmetric  in  k  and  k  . 
This  is  also  the  most  efficient  path  to  produce  H^O  by  at  early  timel.  and  ^ 
since  there  is  a  one  to  one  correspondence  between  [H]  and  [H  O] ,  the 
sensitivity  of  these  concentrations  to  k  and  k„  are  the  same. 

This  system  has  not  been  run  with  a  different  frequency  set  to  deduce 
the  errors  from  interferences.  It  is  estimated,  however,  for  this  set  of  results 


Table  2. 9 


SENSITIVITY  OF  THE  CONCENTRATIONS  AT  10“®  SEC 
TO  UNCERTAINTIES  IN  THE  RATE  COEFFICIENTS 


FOR  THE  H-0  SYSTEM 

Rate  Coefficient 

Sensitivity 

Rate  Coefficient  Sensitivity 

_  _  -  _ 

- (11. 

4 

. 

- - - 

- . 

Insensitive  to 

uncertainties 

Insensitive  to  uncertainties 

in  any  rate  constant. 

in  any  rate  constant. 

_  _  _  _ 

-  -  -  (Hi 

,  /.ll 

xlO  - - 

.  _  _  _ 

- (0]  X  10®® - 

4.67 

>^2 

1.79 

0.202 

-0.016 

0.202 

h 

0.0065 

h 

0.010 

\ 

0.0065 

\ 

0.0023 

2.25x10”® 

h 

0.00165 

>^4 

-1.95x10”® 

8.  OxlO"® 

7.87x10"® 

4.1x10”® 

1.20x10”® 

•'3 

-3. 0xl0“® 

•'7 

-8.46x10“^ 

-  -(OH]  X 

10'® . 

_  _  _  _ 

-  -[HjjO]  X  10®®  -  -  . 

>'8 

5.9 

2.03 

-0.00203 

S 

2.03 

1.76x10”® 

h 

0.121 

•'6 

1.14x10"® 

0.0066 

1.61x10"^ 

>'4 

8.35xl0”‘ 

9.73x10”® 

>^1 

-6.34x10”' 

7.08x10”® 

>'9 

4.22x10' 

-7.43x10”® 

-2,81x10”' 

•'3 

3.13x10”® 

1.33x10”' 

that  no  credence  should  be  given  sensitivities  smaller  than  1  percent  of  the 
maximum  for  that  species. 

2. 5.2. 2  Interpretation  of  Results  for  the  H-0  System  at 
10"3  sec 

-3 

At  t  =  1C  sec  the  system  is  about  to  make  dramatic  changes.  Therefore, 

as  the  rate  constants  are  shifted,  one  would  expect  the  concentrations  to  change 

-3  ' 

dramatically  at  t  =  10  sec.  This  is  indeed  the  case,  and  one  finds  the  various  f 

sensitivity  results  listed  in  Table  2.10.  As  predicted,  reactions  5,  6,  and  7  seem 
to  be  most  important  overall,  with  the  early-time  reactions  1,  2,  and  8  of  con¬ 
siderably  less  importance.  It  is  clear  that  the  controlling  reactions  affect  aU 
the  species  about  equally  vdiich  is  to  be  expected,  since  all  species  appear  in 
those  reactions. 

In  running  the  H-0  system  to  late  times,  difficulties  were  encountered 
with  the  equation  integration  package.  The  timesteps  became  excessively  short, 
and  often  the  input  error  criterion  had  to  be  modified  in  order  to  obtain  a  timely 
solution.  The  difficulties  pointed  up  the  need  to  investigate  more  closely  methods 
for  integrating  the  kinetic  equations  more  economically. 

2. 6  COMPARISON  OF  THE  FOURIER  METHOD  WITH  OTHER  METHODS 
FOR  EVALUATING  THE  MULTIPLE  INTEGRALS 

It  is  possible  to  view  the  Fourier  method  as  a  means  for  approximately 
evaluating  the  multiple  integral  of  an  output  function  over  u  space.  If  incom¬ 
mensurate  frequencies  could  be  used,  the  evaluation  would  be  exact.  By  using 
Int^er  frequencies,  the  number  of  points  required  to  determine  the  Fourier 
amplitude  is  finite,  but  an  error  is  introduced  because  of  the  frequency  inter¬ 
ferences. 

One  way  of  obtaining  the  Fourier  amplitudes  without  the  troublesome  inter¬ 
ferences  would  be  to  evaluate  the  n-dimensional  integral  to  which  the  amplitudes 
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Table  2.10 


SENSITIVITY  OF  THE  CONCENTRATIONS  AT  10~^  SEC 
TO  UNCERTAINTIES  IN  THE  RATE  COEFFICIENTS 
FOR  THE  H-O  SYSTEM 


Kate  Cocmccnt  Sensitivity _ Palo  CocfnCcnt  Sensitivity 


- [Ilj) 

xio® - 

•  [H) 

X  10®  -  -  .  . 
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1.92 

‘‘s 
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“a 
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''a 
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“2 

-o.loa 

“2 

0.061 

“4 
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'‘4 
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0.0176 

''1 

^3 

0.0119 

•'a 

•0.00318 

-  -  -  -  [O3) 

X  10®  -  _ 

lO)  X 

10*® 

•  ^ 

-1.77 

“5 

5.77 

-0.153 

0.420 

“e 

-0.089 

■'i 

0.164 

“a 

-0.070 

■'o 

0.162 

'‘2 

-0.047 

0.145 

“4 

-0.047 

>'2 

0.116 

-0.044 

■'a 

0. 115 

“i 

-0.031 

0.10 

■‘3 

0.0091 

“3 

-0.061 

-  -  ton)  X 

.0^“- . 

(llgOj 

X  10®.  -  -  . 

“a 

4.63 

>'5 

2.52 

“7 

0.291 

0.237 

“a 

0.196 

“0 

0.142 

•‘a 

0,135 

■‘a 

0.114 

“a 

0.116 

0.114 

“2 

0.115 

“2 

0.0718 

0.0964 

“9 

0.0614 

“1 

0. 0876 

“1 

0.0358 

“3 

-0.0359 

“3 

-0.00851 
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correspond  (see  Equation  2. 22)  by  an  alternate  method.  One  method  which  was 
tried  was  based  on  simple  Monte  Carlo  integration.  In  this  section  we  comparq 
the  Fourier  method  and  the  corresponding  Monte  Carlo  calculations.  Another 
equivalent  form  was  investigated,  namely  integral  in  Appendix  C  which  is  an 
n-dimensional  s  integration  equal  to  the  u  integral  transformed  by  the  relations 

Uj  =  sinco.s  ,  i  =  l,  n  (2.39) 

when  the  frequencies  co.  are  those  associated  with  the  Fourier  method.  The 

function  used  in  these  comparisons  was  the  "simulation  function, "  Equation  (3. 3) 
with 


Cj(«o)  =  0.1  ,  C.(0)  =  1 

^i  =i  ,  u°  =  0.5  ,  i  =1,  6 
t  =  1.5 


The  results  for  the  Fourier  sine  amplitudes  of  the  input  frequencies  are 
denoted  A^  ,  i  =  1,  6  .  The  frequencies  were 

(1,  21,  31,  37,  45,  49)  , 

201  points  were  used  on  the  interval  (0,  27r)  (no  aliases  interference  of  order  less 
than  5).  The  amplitudes  were  run  for  5  different  permutations  of  frequencies, 
each  of  which  defined  a  different  path  through  the  space.  The  results  are  shown 
in  Table  2. 11.  The  quantity  in  the  column  labeled  "A"  is  the  average  ampli¬ 
tude,  and  the  "Maximum  Error"  is 

A  -A  I 

_ I  maximum 
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and  serves  as  an  estimate  of  the  errors  introduced  by  interferences. 

The  same  one-dimensional  integral  was  performed  by  Monte  Carlo 
selecting  the  s  points  on  the  interval  (0.  2tt)  at  random,  and  calculating  an 
approximate  integral 

N 

'l  “  2^  XI  “l/i  . 

1=1 

The  results  for  the  amplitudes  for  successive  samples  of  200  points  is  shown  in 
Table  2. 12.  These  results  show  that  the  number  of  points  required  by  this 
Monte  Carlo  procedure  to  give  statistical  accuracy  is  much  greater  than  the 
number  required  in  the  Fourier  analysis.  Stated  another  way,  for  200  points, 

the  degree  of  accuracy  of  the  Fourier  method  far  surpasses  that  of  the  Monte 
Carlo  Integration. 


Another  method  was  to  select  six  values  of  u  at  random  for  each  func¬ 
tion  evaluation,  and  approximate  the  integral  by 


N 


The  results  are  shown  in  Table  2. 13  for  two  successive  sets  of  200  points. 

While  the  error  Is  meeh  less  In  this  case  than  in  performing  the  one-dimensional 
Monte  Carlo  Integral,  the  expeeted  errors  in  a  particular  200  point  sample  still 
greatly  exceeds  the  errors  in  the  Fourier  sampling. 

Monte  Carlo  was  also  used  to  evaluate  the  six-dimensional  integral 
(I3  of  Appendix  C)  in  s  space  choosing  sfac  independent  values  of  s  to  insert 

separately  into  the  integral.  The  results  were  comparable  to  those  in  Table  2. 13 
and  will  not  be  given. 
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3.  OPTIMIZATION  OF  AN  OUTPUT  FUNCTION 
WITH  RESPECT  TO  VARIATIONS  OF  THE 
RATE  COEFFICIENTS 

An  alternate  method  for  investigating  the  variation  of  the  calculated  con¬ 
centrations  with  variations  of  the  rate  constants  was  investigated  in  which  the 
concentrations  were  optimized  over  the  multi-dimensional  domain  of  rate  con¬ 
stants  bounded  by  the  maximum  and  minimum  estimated  uncertainties. 

Several  methods  for  optimizing  functions  of  many  variables  have  been 
weU  automated  and  are  operational  on  the  S^  computer.  A  code  CUK  was 
written  which  solves  the  problem  by  the  method  of  steepest  ascents.  This 
method  proceeds  by  choosing  the  direction  of  the  maximum  gradient  from  a 
given  starting  point,  and  following  that  line  to  an  extremum.  The  process  is 
repeated  by  taking  this  as  a  new  starting  point,  constructing  a  new  line,  following 
it  to  its  extremum,  and  so  forth  until  a  relative  extremum  in  the  many  dimensional 
space  is  fmmd.  A  refinement  to  this  technique  using  a  sophisticated  algorithm  due 
to  Davidon  and  modified  by  Fletcher  and  Powell^^^  and  Stewart is  embodied 
in  subroutines  DMINl  and  DMIN2.  In  the  Davidon  algorithm,  the  search  for  a 
minimum  proceeds  along  a  line  which  does  not  exactly  coincide  with  the  gradient. 

At  each  step  a  new  direction  is  chosen  dependent  not  only  on  the  gradient  at  the 
new  starting  point,  but  also  on  the  previous  history  of  the  search.  For  some 
functions  this  procedure  is  much  more  efficient  than  one  using  the  gradient 
directly.  DMINl  requires  a  subroutine  which  computes  the  gradient,  which  is 
essentially  that  described  below,  see  Section  3.1.  DMIN2  estimates  the  gradient 
by  repeated  evaluation  of  the  function.  Neither  DMINl  nor  DMIN2  are  capable  of 
handling  constraints,  and  an  extension  of  Davidon's  method  to  include  constraints 
has  been  programmed  into  a  subroutine  called  KEELE. 
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3. 1  FORMULATION  OF  THE  GRADIENTS  FOR  THE  METHOD  OF 
STEEPEST  ASCENTS 

Both  DMINl  and  CUK  require  the  gradient  of  the  functions  to  be  extremized. 
For  the  special  class  of  output  functions  consisting  of  species  concentrations  as 
a  function  of  time,  the  gradient  may  be  integrated  along  with  the  rate  equations  to 
obtain  gradients  at  any  particular  time. 


A  set  of  rate  equations  may  be  formally  written 


dC. 
_ 1 

dt 


n 


E 

i=i 


€  .  k. 
1  ] 


n 


(3.1) 


where  the  index  j  label  (separately)  reactions  in  both  forward  and  backward 
directions.  L(j)  is  the  set  of  indices  corresponding  to  the  reactant  concentra¬ 
tions  on  the  destruction  side  of  the  j  A  equation,  is  the  stoichiometric  coeffi¬ 
cient  of  C.  if  C.  is  produced  in  the  reaction,  is  equal  to -1  if  C.  is 
destroyed  by  the  reaction,  and  is  equal  to  0  if  C.  does  not  appear  in  the 
equation  or  appears  with  the  same  stoichiometric  coefficient  on  both  sides  of 
the  reaction.  m(j,jf )  is  the  stoichiometric  coefficient  of  on  the  destruction 
side  of  the  equation. 


The  partial  derivatives  of  the  i^  concentration  with  respect  to  the  p^ 

rate  coefficient  C  are  found  by  differentiating  Equation  (3.1)  with  respect  to  k 
ip  P 


ac, 

“at 


fv  o  L/.  I  I 

2£  =  A  _ L  =  M 

at  ak  i  11 

P  ifL(p) 

n 


(C,) 


m(p,i ) 


m(p*,^ ) 


P  jfcL(p*) 


(3.2) 


i=l  jJ'=l  ’  ^«L(j) 
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p*  labels  the  reverse  of  reaction  p ,  and  the  other  symbols  are  defined  after 

Equation  (3.1).  The  second  term  on  the  ri^t  side  of  Equation  (3.2)  makes  use 

of  the  assumption  that  is  a  constant  over  the  uncertainty  space  (i.e. , 

there  is  no  uncertainty  in  the  equilibrium  constant).  This  set  of  equations  for 

p  =  1 ,  n  forms  a  linear  set  of  equations  in  C.p  which  can  be  integrated  to 

obtain  C  (t)  by  the  same  Gear  algorithm  used  for  integrating  the  basic 

i.P 

rate  equations.  A  version  of  the  KEM  code  has  been  written  to  accomplish  this. 

3.2  TEST  CALCULATIONS 

Two  types  of  test  calculations  were  used  to  investigate  the  various  methods 
of  optimization.  First,  the  properties  of  a  ’’simulation  function”  were  investigated 
to  determine  the  consistency  and  relative  economy  of  the  various  methods. 

Secondly,  the  CUK  method  was  applied  to  the  nitrogen-ojqrgen  (N-0)  chemical 
system. 

•  (4) 

A  ’’simulation  function”  was  suggested  by  Professor  Shuler  which  is  not 
the  solution  to  any  known  set  of  kinetic  equations,  but  which  possesses  some  of 
the  same  characteristics  as  kinetic  solutions.  The  function  is  given  by 

Ci(t)=  |cj(0)  C^(oo)  j^C^(0)+  |cj(«>)  -  C^(0)|  e  ^  j 

\  P  2  2 

Cg(0)C3(0)k^k^t  e 

0^(0)  0^(0) 

2  2 

C4(0)C5(0)+Cg(0)k^kgt  + 

Cj^(t)  is  a  function  of  six  rate  constants  as  well  as  the  time.  The  following  substi¬ 
tution  was  made. 
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o 
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u. 

k.  =  k°  e  ' 
1  1 


(3.4) 


so  that  the  function  is  viewed  as  a  function  of  u  in  the 

i 

.  i  =  1.  6  . 


range 


o  o 

<u,<u, 


The  time  t  was  fixed  at  1. 5  eo  that  only  the  funetional  dependence  in  u 


(3.5) 


remains. 


It  is  desirable  to  confine  the  domain  of  interest  in  u  space  to  be  contained 
withm  some  n-dlmenslonal  box.  Since  some  of  the  methods  are  designed  tor 
unconstrained  systems,  the  above  function  nas  multiplied  by  an  importance  tunc- 
tion  to  form  the  "associated  simulation  function" 


n  2  .  o2 
-u./u, 

1  1 


A' 

P(u)=e  xCj(t) 

where  (t)  is  given  by  Equation  (3. 3). 


(3.6) 


A  comparison  was  made  between  the  four  methods:  CUK.  DMim  DMIN2 

and  KEELE.  Ihe  output  ftmotion  is  the  associated  simulation  function  described  ’ 
above,  with  the  parameters  shown  in  Table  3. 1. 


Table  3. 1 

PARAMETERS  FOR  THE  OPTIMIZATION  METHOD  COMPARISON 

Rate  coefficient  nominal  values:  k.  =  1,  i  =  i,  q 
Rate  coefficient  uncertainties;  u°  =  1 ,  i  =  i,  6 

t  =  1.5 

C.(0)  =1  ,  i  =  1,  e 

Cj(«>)  =0.1 
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The  results  using  the  four  methods  are  shown  in  Table  3. 2;  all  using  comparable 

-5 

minimum  step  sizes  of  10  .  These  results  show  that,  for  this  function  of  six 

variables,  the  methods  appear  to  be  consistent. 

Two  starting  points  were  tried,  and  for  this  output  function,  the  same 
end  point  is  reached.  For  complicated  output  functions,  one  would  not  expect 
the  same  end  point  for  different  starting  points. 

The  cost  of  these  calculations  are  all  proportional  to  both  number  of  func¬ 
tion  evaluations  and  the  number  of  derivative  evaluations.  The  relative  economy 
of  the  CUK  method  over  the  more  sophisticated  methods  led  us  to  concentrate  our 
further  investigation  on  this  method.  For  functions  more  complicated  than  the 
above,  it  is  anticipated  that  the  other  methods  mi^t  well  be  superior. 

The  CUK  code  was  used  on  tlie  associated  simulation  function  for  different 
parameters,  and  with  different  variations  in  the  method.  In  the  sets  of  results 
shown  below,  the  parameters  entering  the  associated  simulation  function  were 
those  given  in  Table  3.3. 


Table  3. 3 

PARAMETERS  ENTERING  THE  ASSOCIATED  SIMULATION 
FUNCTION  IN  TEST  PROBLEMS  BELOW 


k°  =  1.  k°  =  2,  k°  =  3, 
k°  =  4,  k°=5,  k°  =  6 

C.(0)  =  1,  i  =  l,6, 

C^H  =  0.1 


Rate  coefficients  nominal 
values : 


t  =  1.5 


Table  3.2 


RESULTS  OF  COMPARISON  OF  OPTIMIZATION  METHODS 
FOR  THE  ASSOCIATED  SIMULATION  FUNCTION 


Part  A;  Minimum  from  Starting  Point 

u°  =  0. 1.  i  =  1,  6 
1 


CUK 

DMINl 

DMIN2 

KEELE 

Function 

0.  356846 

0.356843 

0. 356755 

0. 356866 

Value 

k. 

1.11918 

1. 11994 

1.12342 

l.i!P39 

1 

k 

0.914977 

0.909268 

0.914298 

0.915617 

2 

k„ 

1.45013 

1.44725 

1.41251 

1.45692 

3 

^4 

0.679556 

0. 669602 

0.670227 

0.677832 

k_ 

1.31851 

1.32621 

1.31683 

1.32631 

5 

k- 

1.12414 

1.12282 

1.12951 

1.12901 

6 

u 

0.112598 

0.113276 

0.116388 

0. 112784 

1 

u 

-0.088856 

-0. 095116 

-0.089598 

-0. 088157 

2 

u 

0.371653 

0. 369666 

0.345370 

0.376326 

3 

"4 

-0.386315 

-0.401071 

-0.400138 

-0.388856 

u 

0.276503 

0.282322 

0.275225 

0.282404 

5 

^6 

0.117014 

0.115841 

0.121783 

0.121343 

No,  of  Function  Calls 

46 

115+ 

184 

49 

No,  of  Gradient  Calls 

10 

25 

24 

53 


0 


Table  3.2  (Contd) 


Part  B:  Minimum  from  Starting  Point 
u°  =-0.1,  i  =  1,  6 


CUK 

DMINI 

CMIN2 

KEELE 

Function 

0. 3  56  845 

Value 

0.356839 

0.356755 

0.356841 

1.12464 

1.13471 

1.12343 

1.12512 

^2 

0.916838 

0.924054 

0.914097 

0.919816 

1.44718 

1.44315 

1.41256 

1.45030 

0.670045 

0.669190 

0.670133 

0.673464 

1.29726 

1.30969 

1.31687 

1.31326 

^6 

1.12978 

1.14590 

1.12904 

1.13609 

^1 

0.117460 

0.126376 

0.116384 

0.117894 

^2 

-0.086825 

-0. 078985 

-0.089819 

-0.083581 

^3 

0.369616 

0.366831 

0.345401 

0.371773 

“4 

0.400410 

-0.401689 

-0.400280 

-0.395321 

0.260252 

0.269792 

0.275258 

0.272509 

0.122021 

0. 136189 

0.121367 

0.127595 

No.  of  Function  Calls 

40 

604 

284 

51 

No.  of  Gradient  Calls 

8 

100 

25 

54 


The  results  for  both  maximization  and  minimization  of  the  associated 

simulation  function  are  shown  in  Table  3.4.  In  Table  3.4,  Part  A,  the  maximum 

uncertainty  in  the  rate  coefficients  is  u.  =  0. 5,  in  Part  B  the  uncertainty  is 

assumed  to  be  u°  =  1 . 

i 

3.2.1  Discussion  of  Table  3.4 

Two  methods  were  used  to  search  for  the  maximum  values  of  the  function, 
constrained  and  unconstrained.  In  the  constrained  method,  values  of  u  were  not 
allowed  to  be  outside  the  box. 


When  a  gfradient  indicated  that  larger  values  of  the  function  occurred  outside  this 
box,  the  direction  of  search  was  taken  along  the  boundary.  The  results  for  the 
constrained  and  ^constrained  maxima  are  therefore  different  for  cases  in  which 
the  maximum  lies  outside  the  box. 

Each  optimum  of  the  function  was  evaluated  for  2  or  3  starting  points 
designated  by  the  row  labeled  initial  point.  The  columns  labeled  ±0. 1  denote  the 
initial  point 

u./u°  =  ±0.1  ,  i  =  l,  6  .  (3.8) 

The  column  denoted  "corner"  represents  an  initial  point  with  coordinates  given  by 
the  footnote  asterisk, .  The  rows  denoted  k.  are  the  6  values  of  the  rate  coeffi¬ 
cients  specifying  the  optimal  values,  and  the  ^re  then  coordinates  in 

u  space  normalized  to  a  unit  hypercube.  In  Part  A,  for  example,  the  maximum 
unconstrained  value  of  the  associated  simulation  function  is  found  at 
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Table  3.4 


RESULTS  OF  OPTIMIZATION  OF  THE 
ASSOCIATED  SIMULATION  FUNCTION  (ASF) 


PartA:  Rate  Constant  Uncertainties 
u°  =  0. 5,  i  =  1,  6 


MINIMUM 

MAXIMUM 

Constrained 

Uneons 

trained 

Initial  Point 

1 

e 

Comer* 

■9 

Corner** 

0.1 

Corner** 

1.336502 

1.336492 

1.336565 

0.942881 

0.942877 

0.942912  . 

0.958178 

0.986023 

“2 

1.994806 

1.994808 

1.9948(K, 

2.038508 

2.038510 

2.038437 

2.038124 

2.54384 

■s 

3.732751 

3.732694 

3.7330,'j8 

1.819592 

1.819592 

1.819592 

1.699020 

2.84773 

“4 

3.867460 

3.867470 

3.867411 

4.777368 

4.777428 

4.777113 

4.720191 

5.716707 

•'8 

5.012121 

5.012122 

5.012120 

4.998960 

4.998960 

4.999181 

4.999231 

2.47980 

“o 

6.000996 

6.000998 

6.000997 

5.999934 

5.999934 

5.999980 

5.999951 

4.99548 

A  0 

“A 

0.58 

0.58 

0.58 

-0.12 

mm 

-0.12 

-0.087 

-0.028 

/  0 
“2/“2 

-0.004 

-0.007 

-0.007 

0.04 

0.04 

0.004 

0.480 

/  0 
“3/“3 

0.44 

0.44 

0.44 

-1.0  - 

-1.0 

-1.13 

-0.1 

-0.067 

-0.069 

-0.070 

0.35 

0.35 

0.36 

0.33 

0.71 

V“”5 

0.005 

0.004 

0.001 

-C.0005 

-0.0004 

-0.0005 

-0.003 

-1.4 

0.001 

-0.001 

-0.001 

-0.00005 

0.00002 

-0.00004 

0.00002 

-0.367 

Optimal  Value 
of  the  ASM 

0.130960 

0. 130968 

0.130968 

0.173318 

0.173318 

0.173318 

0.174365 

0.161049 

No.  of  Function 
Evaluations 

102 

75 

m 

B 

57 

65 

61 

■ 

No.  of  Derivative 
Evahiatlona 

30 

20 

• 

16 

16 

14 

16 

15 

14 

*lnlUal  Point)  u/u°-(0.9,  -0.8,  0.9,  -0.9,  0.8,  -0.8>  **InltlnI  Point)  u/u®-(-0.5,  0.9,  -0.7,  0.95,  -0.95,  -45) 


Table  3.4  (Contd) 


Part  E: 


Rate  Constant 
o 


u . 


=  1.  i 


Uncertainties 

=  1.  6 


MINIMUM 

I  MAXIMUM 

Uncona 

trained 

1  Constrained 

Initial  Point 

IIHEi 

Corner* 

0.1 

Corner* • 

1 

Corner** 

‘‘l 

1.732246 

1.732198 

0.972464 

0.891551 

1  0.984621 

0.885884 

■‘2 

1.988241 

1.988255 

2.204691 

3.334084 

2.228403 

2.814347 

4.172936 

4.172686 

1.103838 

2.881929 

0.873297 

2.93961 

K 

3. 814059 

3.814118 

7.507727 

10.87313 

7.161887 

14.42193 

6. 067699 

5.067709 

4.997028 

1.839397 

4.998524 

1.49980 

“e 

6. 006362 

6.006398 

6.999816 

3.641476 

6.999911 

3.740973 

“i 

0.65 

0.65 

-0.11 

-0.015 

-0.12 

“2 

-0.005 

-0. 005 

0.61 

0.11 

0.34 

“3 

0.33 

0.33 

-0.04 

-1.23 

-0.02 

“4 

-0.047 

-0.047 

0.63 

1.0 

0.68 

1.28 

"6 

0.014 

0.014 

-0.0006 

-1.0 

-0.0004 

-1.2 

”6 

0.002 

0.002 

-0.00002 

-0.49 

-0.0001 

-0.47 

Optimal  Value 
of  the  ASM 

0.126953 

0. 126953 

0.362646 

1.25303 

0.376890 

1.49207 

No.  of  .  -5 

Evalua* 

tiona 

66 

32 

69 

61 

36 

46 

26 

73 

39 

B 

19 

12 

10 

10 

7 

30 

6 

•Initial  Point:  u  -  (0.9i,  -0.95,  0.96,  -0.95,  0.7,  -0.9) 


••Initial  Point:  u  (-0.6,  0.9,  -0.7,  0.95,  -0.96,  -0.6) 


a  point  which  is  outside  the  uncertainty  box,  ^ile  in  the  constrained  case,  u 

_  Q 

lies  on  the  surface  of  the  box,  u  =  -\i  . 

3  3 

In  Part  A,  with  the  estimated  uncertainly  u°  =  0.5,  i  =  1,  6  the  optimum 
function  value  is  found  to  be  independent  of  the  initial  point  of  the  search  for 
both  the  minimum  and  constrained  maximum  cases,  but  the  unconstrained  case 
finds  different  maxima.  In  Part  B,  even  the  maximum  found  in  the  constrained 
search  is  dependent  on  the  initial  point.  This  emphasizes  the  fact  that  such  a 
search  finds  only  a  relative  optimum  \(hich  may  be  path  dependent. 

Finally,  some  words  will  be  said  about  the  numbers  of  function  and  deriva¬ 
tive  evaluations  necessary  to  reach  the  optimum.  These  numbers  depend 
strongly  on  the  complexity  of  the  output  function.  The  method  used  by  the  CUK 
code  is  to  select  a  direction  based  on  a  gradient  calculation,  and  then  search 
along  that  path  with  ever  increasing  distance  until  a  functional  optimum  is 
found  (or,  if  the  search  is  constrained,  until  a  boimdary  is  reached).  At  this 
point,  another  gradient  is  calculated  and  the  search  proceeds  along  the  new  direc¬ 
tion.  The  fact  that  the  number  of  gradient  calculations  is  small  compared  with 
the  number  of  functional  evaluations  indicates  that  the  gradient  varies  relatively 
slowly  in  the  directions  of  search.  Consistent  with  this  observation  are  the  last 

two  rows  in  Table  3.4(B)  which  show  that  increasing  the  minimum  step  size 
-5  -2 

from  10  to  10  (therefore  increasing  the  average  step  size)  reduces  the 
number  of  fimction  evaluations  more  than  it  reduces  the  number  of  gradient  cal¬ 
culations. 

An  interesting  feature  of  the  associated  simulation  fimction  is  that  the 
hyper-corner  in  which  the  maximum  occurs  is  the  reflection  through  the  origin 
of  the  hyper-corner  where  the  minimum  is  found  (there  are  slight  deviations 
from  this  principle  in  the  case  of  u.  in  Part  A  for  both  the  maximum  and 

D 

minimum.  In  each  case,  however,  the  values  of  u„  are  small  so  that  difference 

6 

of  hyper-comer  could  be  the  effect  of  round-off  error. ) 


3.3 


RESULTS  FOR  THE  N-O  CHEMICAL  SYSTEM 


The  N-O  system  investigated  was  the  same  as  that  for  which  test  calcula¬ 
tions  were  run  using  the  Fourier  method.  TTie  set  of  reactions  chosen  was: 


0-+M3220  +  M 

ki. 

N„  +  M5S2N+M 

k-, 

2 

2’ 

2 

NO  +  MiTN  +  O  +  M 

^3’ 

N-  +  O  ^  NO  +  N  k,,  K, 

£>  4  4 

0„  +  N  ^  NO  +  N  k_,  K_ 

A  o  o 


The  forward  rate  constants,  kj ,  and  the  equilibrium  constants  K.  were  chosen 

as  described  in  Table  2.4,  with  the  exception  that  K_  was  set  equal  to  0. 59 

5 

instead  of  59. 0.  The  uncertainty  in  each  rate  constant  was  assumed  to  be  a 
±2 

factor  10  ,  that  is,  the  rate  coefficients  were  assumed  to  lie  in  the  range 


k. 

1 


(3.9) 


-4.606<u.  <4.606  . 

The  concentrations  were  initialized  as  shown  in  Table  3.5,  and  the  system 

-4 

allowed  to  react  isothermally  at  constant  density  until  t  =  10  sec  when  the 
various  species  reach  concentrations  which  were  considered  as  output  functions 
depending  on  the  various  rate  coefficient  uncertainties.  The  CUK  code  was  used 
which  is  strictly  a  steepest  ascents  calculation.  The  hyperplanes  defining  the 
boundaries  according  to  Equation  (3.10)  were  taken  to  absolute  consti’aints  not  to 
be  violated. 

The  results  of  the  optimization  is  summarized  in  Table  3. 5  for  the  maximum 
of  each  concentration,  and  in  Table  3.6  for  the  minimum  of  each  concentration. 


ft 


Table  3. 5 

RESULTS  FOR  THE  MAXIMUM  CONCENTRATIONS 
FOR  THE  N-0  SYSTEM 


o 

Initial 

Concentration 

(moles/cm^) 

Maximum 
Concentration 
(moles/cm^) 
after  10“'*  sec 

Number  of 
Function 
Evaluations 

Number  of 
Gradient 
Evaluations 

Boundary 

States 

o 

02 

2.00x10"® 

1.98  X  10“® 

22 

9 

0 

0 

4.00 X  10“® 

13 

5 

+ - + 

^2 

8.00x10“® 

8.00X  10“® 

22 

9 

1 

N 

0 

8.68xl0“® 

17 

10 

+  -++  + 

G 

NO 

0 

4.70  x  l0“® 

38 

16 

M--M- 

o 

Table  3.6 

RESULTS  FOR  THE  MINIMUM  CONCENTRATIONS 

Q 

FOR  THE  N-O  SYSTEM 

Initial 

Minimum 

Nmnber  of 

Number  of 

Boundary 

Concentration 

Concentration 

Function 

Gradient 

Status 

(moles/cm^) 

(moles/cm^) 

Evaluations 

Evaluations 

0 

after  10"*  sec 

O2 

2.00X  10“® 

2.34  X  10"^^ 

16 

9 

+  +++- 

0 

C 

9.12  X  10"^^ 

17 

9 

o 

^2 

8.00x10"® 

3.69  X  10“® 

25 

10 

+  -++  + 

N 

0 

1.77x10"^® 

14 

7 

NO 

0 

5.61x10"^® 

18 

10 

+ - + 

O 


o 


60 


In  both  tables,  the  columns  need  no  further  explanation  except  for  the 
column  designated  "Boundary  Status."  The  codes  M  have  a  meaning  best 

described  by  example.  The  boundary  status  of  N  in  Table  3.5  is  +-  +  +  +  . 

This  means  k  ,  k.,  and  k  equal  to  their  upper  experimental  uncertainty 
limit  and  k  equals  to  its  lower  experimental  uncertainty  limit  will  produce 
the  maximum  concentration  of  N  at  10  sec.  If  an  M  appears  in  the  i  ft  position 
from  the  left,  the  optimal  '.^lue  of  the  concentration  occurred  for  a  value  of  the 
ith  rate  coefficient  between  its  upper  and  lower  experimental  uncertainties. 

Only  in  the  case  of  the  maximization  of  NO  did  the  extremum  lie  at  an  interior 
point  for  a  particular  rate  constant. 

Another  observation  may  be  made  that  the  hyper -corner  in  which  the 
maximum  is  found  is  not  the  reflection  of  that  for  the  minimum  of  the  same 
species.  In  this  way  this  chemical  system  is  more  complicated  than  the  simula¬ 
tion  function  considered  earlier.  This  may  be  understood  in  terms  of  the 
chemical  kinetics  by  consideration  of  an  example,  the  N  concentration.  In 
Table  3. 6,  the  N  concentration  is  shown  to  be  minimized  by  maximizing  k^  . 
Referring  to  the  reaction  set  on  page  59  this  is  reasonable  since  reaction  5 
contributes  strongly  in  the  forward  direction  since,  by  minimizing  kj  ,  there 
is  an  abundance  of  0_  present.  On  the  other  hand,  Table  3.5  also  shows  that 

A 

the  N  concentrations  will  be  maximized  by  maximizing  k^  .  In  this  case,  kj^ 

is  large  so  that  there  is  an  abundance  of  O  atoms,  and  the  reverse  of  reaction  5 

dominates.  Therefore,  to  increase  k^  increases  the  N  concentration.  It  is 

5 

questionable  if  results  such  as  this  could  have  been  reached  by  chemical  intuition 
alone,  and  we  feel  this  exchange  demonstrates  a  useful  application  of  the  method. 

In  addition  to  the  information  presented  in  Tables  3.5  and  3.6,  one  also 
has  values  of  the  gradients  and  output  function  at  various  points  along  a  path  from 
an  initial  point  and  the  final  optimum  which  can  be  useful  and  give  a  measure  of  the 
relative  sensitivity  of  the  output  function  to  uncertainties  in  the  various  rate 


61 
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o 


o 


o 


o 


') 


coefficients.  The  method  has  the  characteristic  that  the  sampling  of  points 
emphasizes  either  lucge  values  or  small  values  rf  a  particular  output  function 
which  might  be  suitable  for  problems  in  which  either  values  of  a  concentration 
above  or  below  some  threshold  is  of  interest,  for  example,  in  the  investiga¬ 
tion  of  the  onset  of  radar  blackout  in  the  atmosphere  as  a  function  of  various 
rate  coefficients. 

The  optimization  methods  have  the  disadvantage  that  a  separate  analysis 
must  be  made  for  each  output  function  of  interest,  whereas,  in  the  Fourier  method. 

the  sensitivity  of  all  species  in  the  system  may  be  investigated  by  using  the  same 
sample  of  points. 

The  cost  of  the  optimization  procedure  (for  one  concentration)  is  less 
than  the  Fourier  method  in  the  sense  that  fewer  function  evaluations  are  necessary. 
This  effect  is  more  than  balanced,  however,  by  the  increased  cost  required  to 
evaluate  the  gradients  since,  for  n-independent  reaction  rates  the  cost  of  per¬ 
forming  the  time  integration  of  the  derivatives  is  approximately  n^  times  greater 
than  a  calculation  in  which  no  gradients  are  computed.  It  may  be  that  use  of  an 
optimization  method  such  as  DMIN2  or  the  SIMPLEX  method^®)  which  replaces  the 
gradient  calculations  by  differences  of  previously  calculated  output  functions  may 
present  a  more  economical  alternatives 
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4.  SUMMARY  AND  FUTURE  EFFORTS 


O 


Two  types  of  methods  have  been  developed  to  investigate  the  sensitivity 
of  the  results  of  complex  calculations  to  uncertainties  in  the  parameters  entering 
the  calculations,  the  Fourier  method,  and  the  optimization  methods.  For  a  particular 
output  function,  the  Fourier  method  provides  a  sensitivity  number  for  each 
independent  reaction  rate  coefficient  whose  magnitude  characterizes  the  relative 
importance  of  the  uncertainty  in  that  coefficient  in  determining  the  uncertainty  in 
the  output  function.  This  relative  importance  is  an  average  importance  over 
the  many  dimensional  domain  of  the  reaction  rates  consistent  with  the  experimental 
uncertainty  in  each  rate.  The  optimization  methods  investigate  both  the  value  of 
the  output  function  and  its  gradient  with  respect  to  the  rate  coefficients.  The 
latter  is  a  measure  of  local  sensitivity  at  a  point.  The  uncertainty  space  is  not 
sampled  uniformly,  but  regions  with  either  large  or  small  values  of  the  output 
function  are  treated  preferentially.  These  two  methods  have  application  to 
distinctly  different  classes  of  problems. 


For  the  real  chemical  systems  considered  so  far,  there  are  some  features 
of  the  sensitivity  results  that  can  be  predicted  before  the  calculation.  There  are 
others  that  can  not  be  easily  predicted  beforehand,  but  can  be  verified  qualitatively 
by  hindsight.  The  quantitative  results  (which  of  two  important  rate  uncertainties 
is  most  important)  are  ofteh  extremely  hard  to  predict  on  intuitive  grounds  alone, 
and  demonstrate  the  utility  of  the  methods. 

At  present  we  feel  that  both  methods  are  useful  in  answering  different 
questions  about  certain  systems.  Since  the  Fourier  method  is  least  understood, 
and  appears  to  be  most  promising  from  the  standpoint  of  economy,  our  immediate 
plans  are  to  concentrate  on  this  method.  First  priority  is  to  apply  the  method  to  a 
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wider  class  of  non- equilibrium  chemistry  problems,  in  order  to  get  a  feel  for  the 
method.  It  will  be  advantageous  for  large  systems  (depending  on  more  than,  say, 

20  independent  reactions)  to  learn  ways  to  split  the  system  up  into  certain  indepen¬ 
dent  pieces  for  separate  analysis.  Another  interesting  question  \\*ich  will  be 
addressed  is  to  investigate  alternate  ways  to  calculate  the  multiple  integral  form 
of  the  Fourier  coefficient  using  quadrature  rules.  This  would  eliminate  the  trouble¬ 
some  effects  of  interferences  between  frequencies. 


There  are  improved  optimization  techniques  which  do  not  require  the 
evaluation  of  the  gradient,  which  may  make  these  methods  less  expensive  and, 
therefore,  more  useful.  With  these  methods  there  is  always  the  problem  that  an 
optimum  is  a  local  optimum  and  a  global  optimum  requires  many  calculations 
beginning  from  different  initial  points.  The  selection  of  initial  points  is  again  a 
sampling  problem,  and  it  appears  that  considerable  experience  would  be  necessary 
to  adequately  investigate  the  uncertainty  domain. 

Another  area  for  study  is  the  selection  of  integration  routines  which  will 
give  results  more  economicaUy.  This  is  important  since,  for  large  systems, 
the  calculation  of  the  output  function  is  the  major  expense.  Perhaps  accuracy 
could  be  sacrificed  for  speed  in  a  sensitivity  analysis. 

Finally,  the  methods  should  be  applied  to  a  wide  range  of  problems  outside 

the  area  of  non-cquiltbrium  chemistiy,  for  example,  in  the  area  of  social  systems 

modeling,  or  quality  control  modeling  to  determine  their  applicability  to  such 
systems. 
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Appendix  A 


equivalence  of  rational  and 

INTEGER  FREQUENCIES 


Consider 


T 

I  ^sg(sina)  s . sin^.s,  ....  sin  w,  s) 

T-wio  J  J  M 


(A.l) 


where  only  a  ftalle  number,  M .  of  rational  frequenoles  are  present.  By  rational 
frequencies  we  mean  that  they  may  be  represented  by 

<0=^1 

where  p^  and  q^  are  relatively  prime  integers. 

Let 


M 


q  = 


j=i 


and 


qwj  =  a,j  , 


then  quj-iuj  is  an  integer,  and  Equation  (A.l)  may  be 


rewritten 


lim 

T-*«o 


I 


ds  g  (sin<o^ s/q,  ,  sin  s/q)  . 


(A.  2) 


But  since  g  is  a  periodic  function  of  s'  =  s/q  with  period 


2tt  since  all  the 
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wj  are  integers,  this  gives 


1  = 


f2n 


ds' g(sin  w!  s’,  ....  sin  cu’  s')  . 


The  infinite  integral  in  Equation  (A.l)  with  rational  frequencies 
to  a  finite  integral  with  integer  frequencies. 


(A.  3) 

is  therefore  equal 
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Appendix  B 


EXTRACTION  OF  FOURIER  AMPLITUDES 
BY  A  FINITE  SUM 


The  purpose  of  this  Appendix  is  to  show  that  one  can  evaluate  the  Fourier 
amplitudes  of  a  periodic  function  f(s)  by  the  following  prescription 

(^1  2 

I  I  N  2^  t  1  ^  amplitudes  (B.  1) 

'  ®n  M 

where  A^  is  defined  in  Equation  (2, 7),  and  B^  is  the  corresponding  cosine  ampli¬ 
tude.  n  is  an  integer  >0  and  f(2TT)  =  f(0).  We  shall  provr^is  relation  for 
A^,  the  sine  amplitudes  only,  although  the  proof  can  be  directly  extended  to  the 
cosine  amplitudes,  B^^ .  Consider 


(B.2) 


where 


<10 

f(s)  =  B  +  (A  sin  ms  +  B  cos  ms)  . 

o  m  m 


Substituting  f  into  A*  one  must  evaluate  a  sum  of  terms 

n 


,  m  '  I  22  A„  sto(2W/N)  . 


(B.3) 
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N 


A* 
n,  m 


j?=i 


f  2Tri(n+ni)i/N  2Tri(n-m)i/N 


-  e 


m 


(B.4) 


-  Re 
N 


^2Tri(n+m)  ^ 


2m(n-in) 


^27ii(n+m)/N  ^ 


2Tii(n-m)/N 


-1 


m 


where  we  have  used  the  formula 


N 


^  r^  =  (r^-  l)/(r-l) 


(B.5) 


i=l 


In  the  brackets  each  term  vanishes  unless  the  denominator  vanishes  since 
2m 


=  1  for  p  integer.  Only  special  values  of  m  relative  to  n  will  there- 

antribute  to  A*, 
n 

bracket,  respectively, 


fore  contribute  to  A* ,  and  those  are,  for  the  first  and  second  terms  in  the 


and 


m  +  n  =  qN  condition  1  , 
|m-n|  =  qN  condition  2  , 


(B.6) 

(B.7) 


where  q  is  an  integer  >  0 .  For  q  =  0,  only  m  =  n  contributes.  For  q  >  1 , 

let  the  special  values  of  m  satisfying  conditions  1  and  2  be  denoted  m|(q,  n) 

and  m*(q,n),  respectively.  Then,  for  N>n,  m*(q,n)  will  be  distinct  from 

n)  (n  =  0  is  an  expection  which  is  uninteresting),  and  every  such  m* ,  and 

m*  will  contribute  to  A*.  Therefore 
2  n 

\  '^n  mj(q,  n)  ^m*(q,  n)]  * 

q=l 
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The  sum  in  the  above  expression  may  be  viewed  as  an  error  term  to  the  evalua¬ 
tion  .  The  amplitudes  in  the  error  term  are  called  "alias  amplitudes"  since 
an  amplitude  for  frequency  m*(q,n)  masquarades  as  one  for  frequency  n  .  For 
the  error  terms  become  small. 

As  discussed  in  Section  2.3.2  the  input  frequencies  are  chosen  to  be 
odd  integers  which  have  no  interactions  of  order  <5  which  coincide  with 
the  input  frequencies.  Since  the  alias  frequencies  are  equivalent  to  other 
frequencies,  it  is  consistent  to  choose  N  sufficiently  great  that  there  are  no 
aliases  of  interferences  of  order  <5  which  coincide  with  input  frequencies  either. 
Thus,  N  is  chosen  by  trial  and  error  so  that  condition  1  is  satisfied, 

"i  +  t^m  ^  qN  (B.  9) 

where  03.  is  an  input  frequency  and  03  are  all  input  frequencies  and  inter- 
1  m 

ferences  of  order  less  than  or  equal  to  4.  One  can  show  that  a  sufficient 
condition  to  satisfy  Equations  (B.  6)  and  (B.  7)  is 


N  >4o; 


max 


N 

even 


(B.  10) 


but  by  trial-and-error  values  of  N  slightly  smaller  than 


4  cu  may  be  found, 
max 


70 


Appendix  C 


A  THEOREM 

In  this  Appendix  we  present  a  theorem  which  relates  the  Fourier  amplitude 
of  a  function  to  an  average  of  the  output  function  over  the  n-dimensional  space 
spanned  by  the  rate  coefficient  uncertainties. 

Theorem;  Assume  f(u  ,  u  ,  ....  u  )  is  a  polynomial  of  order  p-1  in  the 

la  n 

u's,  over  the  domain 


o  o 

<«!<«! 


o  o 

-u  ^  u  <  u 
n  n  -  n 


Assume  also  a  set  of  integer  input  frequencies  w, ,  i  =  1 ,  \(iiich  are  chosen 
so  that  no  interferences  of  the  order  of  less  than  p  coincide  with  any  input 
frequency.  Then  the  Fourier  amplitude  corresponding  to  the  Lto  input  frequency 
satisfies  the  following  relation: 


■if 


ds  sin  0)^8  f(u^(s),  u^Cs)  ...  u^(s)) 


u  u 

^  1  • 

^  f 

Ti  n  r 


f”  ^  “l  -  V 

*  *  *  I  n  f  ' 


(C.l) 
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yiiere 


Proof: 


/V  o  • 

Uj(s)  =  sin  Wj  s  . 


For  convenience  the  two  integrals  in  Equation  (C.  1)  will  be  referred  to  as 
and  I  ,  respectively.  The  first  step  is  to  transform  I  by  the  transformation 

M  2 


u,  =  u,  sin  cj,  s. 
i  i  i  1 


i  =  1,  n 


(C.2) 


then 


2;r  2;r  2ir 

=  ““if  . J  ^“n"‘”VL'<VV-”W 

A  A  A 


si 


3  • 


The  next  step  is  to  show  that  I-  -  I 

o  X 


Since  f  is  a  polynomial  of  degree  p-1 ,  we  can  expand  it  in  a  Maclaurin 

series 

. “n>  "o  “i  ^  X! 

i  i.j 

where  the  last  term  contains  a  product  of  p-1  u's .  Substituting  Equation  (C.2) 
for  the  u^ ,  i  =  1,  n ,  one  obtains  a  series  containing  products  of  the  form 

”l  "2 

sin  (Wj^s^lsin  (WgSg)  ....  (C.5) 

These  terms  may  be  expanded  using  the  relation 

Performing  the  s^  integrations  one  t!:.us  obtains 


(C.3) 


'3=“L(f  4E'L„“i'4E  <' 

i=l  i.j=l  / 

where  all  terms  of  f  which  are  differentiated  an  odd  number  of  time  in 

and  an  even  number  of  times  in  all  other  u's  contribute  through  derivatives 
of  order  p-1 . 


If  one  makes  the  substitution 


u.  =  u,  sm  (t).  8 
1  1  1 


(C.7) 


into  the  Equation  (C,4)  and  substitutes  that  into  I^,  one  finds  terms  of  the  form 


sin  (w^s)sin(w2  s) 


(C.  8) 


which  is  slightly  different  from  Equation  (C.5)  in  that  the  same  s  appears  in 
each  term  in  the  product.  The  res’ilt  is  that  when  the  expansion  is  made  to 
reduce  such  a  terin  to  one  linear  in  the  sines  or  cosines,  frequencies  appear 
which  are  interferences  of  the  input  frequencies.  In  fact,  since  f  is  of 
order  p-1  in  the  u's,  interference  frequencies  which  are  of  order  p-1 
will  appear  in  the  spectrum  of  f .  But  the  frequency  set  was  hypothesized  to 
be  free  of  p-1  to  order  interferences  coinciding  with  an  input  frequency  such  as 
Therefore,  only  toe  terms  in  f  which  are  of  odd  order  in  u^^  and  even 
order  in  aU  other  u’s  will  contribute.  But  that  is  the  same  set  of  terms  as 
appears  in  Ig  .  Therefore, 


^1  ~  ^3  • 


(C.9) 


In  the  case  that  f  is  of  higher  order  than  p-1  in  the  u’s,  one  can 
estimate  the  error  between  and  .  (I^  =  Ig  independent  of  frequencies. )  For 
definiteness,  choose  p  =  5 ,  viiich  characterize  the  frequency  sets  used  in  the 
applications  in  Section  2,  and  assuine  f  has  a  fifth  order  component 


f  »  U.IL  u«u  u 
jkJfimq  j  if  m  q 


(C.  10) 


with  j,  k,  i ,  m,  q  all  different.  Substituting  according  to  Equation  (C.  2)  into  f 
in  Ig  gives  zero  contribution  to  Ig  .  On  the  other  hand,  if  one  substitutes  for 
this  fifth  order  term  into  according  to  Equation  (C.  7),  one  can  expand  the 
product  (C.  10)  to  be 


^jkimq 


16 


sin  (Wj  +  cog  +Wg  +  "4  +  W5)s  -  sin(Wj^  ^2  “s  "4 
+  sin  (w^  +  ^^2  "3  ■  "4  ‘  "5^®  ■  ®^  ^"1  "2  "3  “"4  "5^® 

+  ...1 


(C.11) 


If  there  exist  frequencies  w.,  w, ,  . . . ,  oj  such  that  one  of  the  16  interference 

J  Ic  Q 

frequencies  in  Equation  (C.  11)  coincides  with  »  ttie  contribution  from  the  fifth 

Jj 

order  term  will  contribute  to  1^^ .  Since  the  order  of  this  contribution  is 


O  (f^Vs  ^^^)  , 


we  can  write  for  a  general  function  f , 


II  =  I2  +  0(f^V2^"^^ 


(C.12) 
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